00001 #include "epsolver.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 "nssolver.h"
00009 #include "springel.h"
00010 #include "intpoints.h"
00011 #include "vector.h"
00012 #include "matrix.h"
00013 #include "mathem.h"
00014 #include "gmatrix.h"
00015 #include "tensor.h"
00016 #include "vecttens.h"
00017 #include <math.h>
00018 #include <string.h>
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 void solve_epressure ()
00030 {
00031 long i;
00032 double *rhs, *lhs;
00033
00034
00035 for (i=0;i<Lsrs->nlc;i++){
00036
00037 rhs=Lsrs->give_rhs (i);
00038 lhs=Lsrs->give_lhs (i);
00039
00040 nullv (rhs+i*2*Ndofm, 2*Ndofm);
00041 if (Mp->tepsol < epplast_sol)
00042 {
00043 Mb->lc[i].assemble (i,rhs+i*Ndofm,NULL,1.0);
00044 Mb->dlc[i].assemble (rhs+i*Ndofm, NULL,lhs+i*Ndofm);
00045 }
00046 else
00047 {
00048 Mb->lc[i*2+0].assemble (i*2+0,rhs+(i*2+0)*Ndofm,NULL,1.0);
00049 Mb->lc[i*2+1].assemble (i*2+1,rhs+(i*2+1)*Ndofm,NULL,1.0);
00050 }
00051
00052 epressure_solver (i);
00053 }
00054 }
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 void epressure_solver (long lcid)
00068 {
00069 switch (Mp->tepsol)
00070 {
00071 case gep_sol:
00072 general_epressure (lcid);
00073 break;
00074 case gepvarsup_sol:
00075 general_epressure_varsup (lcid);
00076 break;
00077 case epplast_sol:
00078 femplast_epressure (lcid);
00079 break;
00080 default:
00081 print_err("unknown solver of nonlinear equation system is required", __FILE__, __LINE__, __func__);
00082 }
00083 }
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096 void general_epressure (long lcid)
00097 {
00098 long i, j, n, ni;
00099 double zero, err, norf, dl;
00100 double *r,*f,*fi,*fb, *fc, lambda, rlam;
00101
00102
00103
00104 n = Ndofm;
00105
00106 ni = Mp->niep;
00107
00108 zero = Mp->zero;
00109
00110 err = Mp->errep;
00111
00112 dl = Mp->stepep;
00113
00114 fb = new double [n];
00115 fi = new double [n];
00116 memset (fb,0,n*sizeof(double));
00117 memset (fi,0,n*sizeof(double));
00118
00119
00120 r = Lsrs->give_lhs (lcid);
00121 f = Lsrs->give_rhs (2*lcid);
00122 fc = Lsrs->give_rhs (2*lcid+1);
00123
00124
00125
00126
00127 lambda = 0;
00128 for (i=0;i<ni;i++){
00129
00130 print_init(i, "wt");
00131
00132 nullv (f+lcid*n, n);
00133 Mb->lc [lcid].assemble (0,f+2*lcid*n,NULL,1.0);
00134 Mb->dlc[lcid].assemble (f+2*lcid*n,r+lcid*n, NULL);
00135
00136 for (j=0;j<n;j++){
00137 f[j] = (f[j] - fb[j]);
00138 if (lambda < 1.0)
00139 f[j] *= dl/(1.0-lambda);
00140 fc[j] = fb[j];
00141 fb[j] += f[j];
00142 }
00143 lambda += dl;
00144 if (lambda > 1.0)
00145 lambda = 1.0;
00146
00147 rlam = arclengthrv(lcid, 1.0, Mp->rlerrep);
00148 if (fabs(rlam-1.0) > Mp->rlerrep)
00149 {
00150 print_err("arclength could not reach required lambda", __FILE__, __LINE__, __func__);
00151 abort();
00152 }
00153
00154
00155 print_step(lcid, i+1, lambda, fb);
00156
00157 nullv (f+lcid*n, n);
00158 Mb->lc [lcid].assemble (0,f+2*lcid*n,NULL,1.0);
00159 Mb->dlc[lcid].assemble (f+2*lcid*n,r+lcid*n, NULL);
00160 for (j = 0; j < n; j++)
00161 fi[j] = f[j] - fb[j];
00162 norf = ss(fi, fi, n);
00163 norf /= ss(f, f, n);
00164 fprintf (stdout, "\n\n MAIN ITERATION LOOP: NORF=%e, lambda = %e\n", norf, lambda);
00165 fprintf (stdout, "********************************************\n\n");
00166 if ((norf < err) && (lambda == 1.0))
00167 break;
00168 print_close();
00169 }
00170 print_close();
00171 delete [] fi; delete [] fb;
00172 }
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185 void general_epressure_varsup (long lcid)
00186 {
00187 long i, j, n, ni;
00188 double zero, err, norf;
00189 double *r,*rb,*f,*fi,*fb, lambda;
00190 double **br;
00191
00192
00193
00194 n = Ndofm;
00195
00196 ni = Mp->niep;
00197
00198 zero = Mp->zero;
00199
00200 err = Mp->errep;
00201
00202 rb = new double [n];
00203 fb = new double [n];
00204 fi = new double [n];
00205 br = new double*[Mt->nn];
00206 memset (rb,0,n*sizeof(double));
00207 memset (fb,0,n*sizeof(double));
00208 memset (fi,0,n*sizeof(double));
00209 memset (br,0,Mt->nn*sizeof(*br));
00210 for (i = 0; i < Mt->nn; i++)
00211 br[i] = new double[Gtm->give_ndofn(i)];
00212
00213
00214 r = Lsrs->give_lhs (lcid);
00215 f = Lsrs->give_rhs (lcid);
00216
00217
00218
00219
00220 for (i=0;i<ni;i++)
00221 {
00222 print_init(i, "wt");
00223
00224 nullv (f+lcid*n, n);
00225 Mb->lc [lcid].assemble (0,f+lcid*n,NULL,1.0);
00226 Mb->dlc[lcid].assemble (f+lcid*n, NULL,r+lcid*n);
00227
00228 for (j=0;j<n;j++){
00229 fb[j] = f[j];
00230 }
00231
00232 memset (r,0,n*sizeof(double));
00233 lambda = arclengthrv(lcid, 1.0, Mp->rlerrep);
00234
00235
00236 internal_forces (lcid,fi);
00237 for (j = 0; j < n; j++)
00238 f[j] = fb[j] - fi[j];
00239 print_step(lcid, i+1, lambda, fb);
00240 norf = ss(f, f, n);
00241 fprintf (stdout, "\n MAIN ITERATION LOOP: NORF=%e", norf);
00242 if ((norf < err) && (i >= Mp->nselnodep))
00243 break;
00244 if (i < Mp->nselnodep)
00245 {
00246
00247 Gtm->restore_codnum();
00248 fprintf(stdout, "\n removing support in the node number %ld", Mp->selnodep[i]+1);
00249 for (j = 0; j <= i; j++)
00250 {
00251
00252 Gtm->gnodes[Mp->selnodep[j]].cn[0] = 1;
00253 }
00254
00255 n = Ndofm = Gtm->gencodnum();
00256
00257 for (j = 0; j < Mt->nn; j++)
00258 noddispl (lcid, br[j], j);
00259
00260
00261 delete Lsrs;
00262 delete Smat;
00263 Smat = NULL;
00264 delete [] fi;
00265 delete [] fb;
00266 delete [] rb;
00267 Lsrs = new lhsrhs;
00268 Lsrs->alloc ();
00269 rb = new double [n];
00270 fb = new double [n];
00271 fi = new double [n];
00272 memset (rb,0,n*sizeof(double));
00273 memset (fb,0,n*sizeof(double));
00274 memset (fi,0,n*sizeof(double));
00275
00276 r = Lsrs->give_lhs (lcid);
00277 f = Lsrs->give_rhs (lcid);
00278 restore_displ(lcid, br);
00279 }
00280 print_close();
00281 }
00282 print_close();
00283 delete [] fi; delete [] fb; delete [] rb;
00284 for (i = 0; i < Gtm->nn; i++)
00285 delete [] br;
00286 delete [] br;
00287 }
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 void femplast_epressure (long lcid)
00308 {
00309 long i,j,idxp, idxs,ni, *nod, ndofn, *cnum;
00310 double *rhs, *rb, *react, *fpb, *fsb, *forig, lambda;
00311 vector ifor;
00312 ivector cn;
00313
00314 vector sig;
00315 matrix sigt(3,3);
00316 double *p = new double[Mm->tnip];
00317 double *porig = new double[Mm->tnip];
00318 char fn[1001];
00319 char fn2[1001];
00320 FILE *femcad;
00321
00322 ni = Mp->niep+1;
00323
00324
00325
00326 rb = new double[Ndofm];
00327 fpb = new double[Ndofm];
00328 fsb = new double[Ndofm];
00329 forig = new double[Ndofm];
00330 nod = new long [Mp->niep];
00331 react = new double[Mp->niep];
00332 cnum = new long [Mp->niep];
00333 for (i=0;i<ni;i++)
00334 {
00335 print_init(i+1, "wt");
00336 if (i < 1)
00337 fprintf (stdout, "\n MAIN ITERATION LOOP: step %ld\n", i+1);
00338 else
00339 fprintf (stdout, "\n MAIN ITERATION LOOP: step %ld, elem. %ld was removed\n", i+1, Mt->ne+1);
00340 fprintf(stdout, " ---------------------------------------\n\n");
00341 rhs=Lsrs->give_rhs (lcid);
00342 idxs = (lcid*2+1)*Ndofm;
00343 idxp = (lcid*2+0)*Ndofm;
00344 if (i == 0)
00345 {
00346 copyv (rhs+idxp, fpb, Ndofm);
00347 copyv (rhs+idxs, fsb, Ndofm);
00348 addv(rhs+idxp, fsb, Ndofm);
00349 nullv (rhs+idxs, Ndofm);
00350 nullv (forig, Ndofm);
00351 }
00352 if (i > 0)
00353 {
00354 nullv (rhs+idxp, Ndofm);
00355 for (j = 0; j < Ndofm; j++)
00356 {
00357 if (i == 1)
00358 rhs[idxp+j] = (Mp->loadcoef[i-1]-1.0)*fsb[j];
00359 else
00360 rhs[idxp+j] = (Mp->loadcoef[i-1]-Mp->loadcoef[i-2])*fsb[j];
00361 }
00362
00363
00364 }
00365 lambda = newton_raphsonep(lcid, forig);
00366
00367 for (j=0; j < Ndofm; j++)
00368 forig[j] += lambda * (rhs+idxp)[j];
00369
00370 Mm->computenlstresses(0);
00371 for (j=0; j<Mm->tnip; j++)
00372 {
00373 allocv(Mm->ip[j].ncompstr, sig);
00374 Mm->givestress (0,j,sig);
00375 vector_tensor (sig, sigt, Mm->ip[j].ssst, stress);
00376 p[j] = first_invar (sigt)/3;
00377 destrv(sig);
00378 }
00379 if (i == 0)
00380 copyv(p, porig,Ndofm);
00381 else
00382 {
00383 for (j=0; j<Mm->tnip; j++)
00384 {
00385 p[j] -= porig[j];
00386 }
00387 sprintf(fn,"%s%s.jk0", Mp->path, Mp->filename);
00388 sprintf(fn2,"%s%s.jk0.%ld", Mp->path, Mp->filename,i);
00389 femcad = fopen(fn, "at");
00390 char dlcid[50];
00391 sprintf(dlcid, "%ld", i+1);
00392 write_elemscalar(femcad, p, "p", dlcid);
00393 fclose(femcad);
00394 rename (fn, fn2);
00395 }
00396
00397
00398 if (i == Mp->niep)
00399 break;
00400 nod[i] = Gtm->gelements[Gtm->ne-1].nodes[0];
00401 ndofn = Mt->give_ndofn(nod[i]);
00402 allocv(ndofn, ifor);
00403 allocv(ndofn, cn);
00404 Spring->res_internal_forces (lcid, Mt->ne-1, ifor);
00405 react[i] = ifor[0];
00406 Mt->give_node_code_numbers (nod[i],cn.a);
00407 cnum[i] = cn[0];
00408
00409 Mt->ne--;
00410 Gtm->ne--;
00411 destrv(ifor);
00412 destrv(cn);
00413 print_close();
00414 }
00415 print_close();
00416 delete [] rb;
00417 delete [] fpb;
00418 delete [] fsb;
00419 delete [] nod;
00420 delete [] react;
00421 delete [] cnum;
00422 }
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436 void restore_displ(long lcid, double **bckr)
00437 {
00438 long i, j, ii;
00439
00440 for (i = 0; i < Gtm->nn; i++)
00441 {
00442 for (j = 0; j < Gtm->gnodes[i].ndofn; j++)
00443 {
00444 ii=Gtm->gnodes[i].cn[j];
00445 if (ii > 0) Lsrs->lhs[lcid*Ndofm+ii-1] = bckr[i][j];
00446 }
00447 }
00448 }
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 double arclengthrv (long lcid, double rlambda, double rlerr)
00466 {
00467 long i,j,k,n,ni,ini,stop,modif,li, numr;
00468 double a0,a1,a2,l1,l2,dl,dlmax,dlmin,psi,lambda,blambda,dlambda,ddlambda,norf,norfp,norv,norfa,zero,ierr;
00469 double lambdao,ss1,ss2,ss3,ss4,ss5;
00470 double *r,*ra,*ddr,*u,*v,*f,*fa,*fp,*fc,*fi;
00471 long cn, ecn[6];
00472 char file[255];
00473 const char *mode;
00474
00475 FILE *grout;
00476
00477
00478
00479
00480 n = Ndofm;
00481
00482 ni = Mp->nlman->nial;
00483
00484 ini = Mp->nlman->niilal;
00485
00486 zero = Mp->zero;
00487
00488 ierr = Mp->nlman->erral;
00489
00490 dl = Mp->nlman->dlal;
00491
00492 dlmax = Mp->nlman->dlmaxal;
00493
00494 dlmin = Mp->nlman->dlminal;
00495
00496 psi = Mp->nlman->psial;
00497
00498
00499 r = new double [n];
00500 ddr = new double [n];
00501 u = new double [n];
00502 v = new double [n];
00503 f = new double [n];
00504 fa = new double [n];
00505 fi = new double [n];
00506
00507
00508 ra = Lsrs->give_lhs (lcid);
00509 fp = Lsrs->give_rhs (lcid*2);
00510 fc = Lsrs->give_rhs (lcid*2+1);
00511
00512 if (Mp->nlman->hdbackupal==2){
00513 arclopen (li,n,lambda,dl,ra,fp);
00514 }
00515 else{
00516 lambda=0.0;
00517 lambdao=0.0;
00518 li=0;
00519 }
00520
00521 if (li == 0)
00522 mode = "wt";
00523 else
00524 mode = "at";
00525
00526 sprintf(file, "%s%s.epg", Mp->path, Mp->filename);
00527 grout = fopen (file, "wt");
00528
00529 for (j=0;j<n;j++)
00530 fa[j]=fc[j]+(lambda)*fp[j];
00531
00532
00533 norfp = loadincr(fp, lambda, n, Mp->nlman);
00534 modif=0;
00535
00536
00537 long kk;
00538 matrix sm(6,6);
00539
00540
00541
00542
00543 for (i=li;i<ni;i++) {
00544
00545 fprintf(grout, "Lambda %e\n", lambda);
00546 fprintf(grout, "Posuny ux --------\n");
00547 for (kk = 0; kk < Mt->nn; kk++)
00548 {
00549 cn = Mt->give_dof(kk,0)-1;
00550 fprintf(grout, " %2ld %e\n", kk+1, ra[cn]);
00551 }
00552 fprintf(grout, "Zatizeni fx --------\n");
00553 for (kk = 0; kk < Mt->nn; kk++)
00554 {
00555 cn = Mt->give_dof(kk,0)-1;
00556 double fpv = fp[cn];
00557 fprintf(grout, " %2ld %le\n", kk+1, fpv);
00558 }
00559 fprintf(grout, "Tuhosti a sily--------\n");
00560 for (kk = 0; kk < Mt->ne; kk++)
00561 {
00562 if (Mt->give_elem_type(kk) == spring_1)
00563 {
00564 Spring->res_stiffness_matrix (kk,sm);
00565 Mt->give_code_numbers (kk, ecn);
00566 if (ecn[0] > 0)
00567 fprintf(grout, " %2ld %e %e\n", kk+1, sm[0][0], sm[0][0]*ra[ecn[0]-1]);
00568 }
00569 }
00570 fprintf(grout, "\n");
00571 fflush(grout);
00572
00573 copyv (ra, r, n);
00574
00575 blambda=lambda;
00576
00577 fprintf (stdout,"\n arc-length: increment %ld lambda %e dl %e",i,lambda,dl);
00578
00579
00580
00581
00582
00583
00584 stiffness_matrix (lcid);
00585
00586
00587 copyv (fp, f, n);
00588
00589
00590 Mp->ssle->solve_system (Gtm,Smat,v,f,Out);
00591
00592 norv = displincr (v,v,n,lcid);
00593
00594 dlambda = dl/sqrt(norv+psi*psi*norfp*norfp);
00595
00596 for (j=0;j<n;j++){
00597 ddr[j]=dlambda*v[j];
00598 ra[j]+=ddr[j];
00599 fa[j]=fc[j]+(lambda+dlambda)*fp[j];
00600 }
00601 ddlambda=dlambda;
00602
00603
00604 internal_forces (lcid,fi);
00605 subv(fa, fi, f, n);
00606 norf=normv(f,n);
00607 norfa = normv(fa, n);
00608 norf /= norfa;
00609
00610 if (Mespr==1) fprintf (stdout,"\n %e %e norf %e",lambda,dl,norf);
00611
00612 if (norf<ierr){
00613
00614
00615
00616 modif++;
00617
00618 lambda+=dlambda;
00619 if (((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
00620 {
00621
00622 dl/=2.0;
00623 if (dl<dlmin)
00624 {
00625 dl=dlmin;
00626 break;
00627 }
00628 modif = 0;
00629 if (Mespr==1)
00630 {
00631 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf);
00632
00633 }
00634
00635 for (j=0;j<n;j++)
00636 ra[j]=r[j];
00637
00638 lambda=blambda;
00639 }
00640 if (modif>1) {
00641
00642 dl*=2.0;
00643 if (dl>dlmax)
00644 dl=dlmax;
00645 modif=0;
00646 if (Mespr==1){
00647 fprintf (stdout,"\n arc length must be modified (dl increase) dl=%e norf=%e",dl,norf);
00648
00649 }
00650 }
00651 if (fabs(lambda-rlambda) <= rlerr)
00652 break;
00653 Mm->updateipval ();
00654
00655
00656 print_step(lcid, i, lambda, fa);
00657 print_flush();
00658 }
00659 else{
00660
00661
00662
00663 stop=0;
00664 for (j=0;j<ini;j++){
00665
00666
00667 Mp->ssle->solve_system (Gtm,Smat,u,f,Out);
00668 copyv (ddr, f, n);
00669 addv(ddr, u, n);
00670
00671 quadeqcoeff (ddr,v,n,ddlambda,psi,norfp,dl,a0,a1,a2,lcid);
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689 numr = solv_polynom_2(a2, a1, a0, l1, l2);
00690 switch (numr)
00691 {
00692 case -1:
00693 print_err("infinite number of solution of constrained condition in function arclength", __FILE__, __LINE__, __func__);
00694 break;
00695 case 0:
00696 print_err("nonsolvable constrained condition in function arclength", __FILE__, __LINE__, __func__);
00697 break;
00698 case 1:
00699 dlambda = l1;
00700 break;
00701 default:
00702 break;
00703 }
00704 ss1=0.0; ss2=0.0;
00705 ss3=0.0; ss4=0.0; ss5=0.0;
00706 for (k=0;k<n;k++){
00707 ss1+=(ddr[k]+l1*v[k])*f[k];
00708 ss2+=(ddr[k]+l2*v[k])*f[k];
00709 ss3+=(ddr[k]+l1*v[k])*(ddr[k]+l1*v[k]);
00710 ss4+=(ddr[k]+l2*v[k])*(ddr[k]+l2*v[k]);
00711 ss5+=f[k]*f[k];
00712 }
00713
00714 if (ss1/ss3/ss5>ss2/ss4/ss5)
00715 dlambda=l1;
00716 else
00717 dlambda=l2;
00718
00719 for (k=0;k<n;k++){
00720 ddr[k]+=dlambda*v[k];
00721 ra[k]+=u[k]+dlambda*v[k];
00722 fa[k]+=dlambda*fp[k];
00723 }
00724 ddlambda+=dlambda;
00725
00726
00727
00728 internal_forces (lcid,fi);
00729 subv(fa, fi, f, n);
00730 norf=normv(f,n);
00731 norfa = normv(fa, n);
00732 norf /= norfa;
00733
00734 if (Mespr==1){
00735 fprintf (stdout,"\n increment %ld inner loop %ld norf=%e",i,j,norf);
00736
00737 }
00738 if (norf<ierr){
00739 lambda+=ddlambda;
00740
00741 print_step(lcid, i, lambda, fa);
00742 print_flush();
00743 stop=1;
00744 if (((lambda - rlambda) > 0.0) && (fabs(lambda-rlambda) > rlerr))
00745 stop=0;
00746 if (fabs(lambda-rlambda) <= rlerr)
00747 stop=2;
00748 if (stop > 0)
00749 Mm->updateipval();
00750 break;
00751 }
00752 }
00753
00754
00755 if (stop == 2)
00756 break;
00757
00758 modif=0;
00759 if (stop==0){
00760
00761 dl/=2.0;
00762 if (dl<dlmin){
00763 dl=dlmin;
00764 break;
00765 }
00766 if (Mespr==1){
00767 fprintf (stdout,"\n arc length must be modified (dl decrease) dl=%e norf=%e",dl,norf);
00768
00769 }
00770
00771 copyv (r, ra, n);
00772
00773 lambda=blambda;
00774 }
00775 }
00776 if (stop==0)
00777 continue;
00778
00779
00780 if (i%10 && Mm->plast)
00781 continue;
00782
00783
00784 }
00785
00786
00787
00788
00789
00790 if (Mp->nlman->hdbackupal==1)
00791 arclsave (i,n,lambda,dl,ra,fp);
00792
00793
00794
00795
00796
00797 delete [] r;
00798 delete [] fi;
00799 delete [] fa;
00800 delete [] f;
00801 delete [] v;
00802 delete [] u;
00803 delete [] ddr;
00804
00805
00806 return (lambda);
00807 }
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
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
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198 double newton_raphsonep (long lcid, double *pfp)
01199 {
01200 long i,j,k,n,ni,ini;
01201 double lambda,blambda,dlambda,dlambdamin,zero,err,norf,norfa;
01202 double *r,*rb,*dr,*f,*fi,*fb,*fc,*fp;
01203
01204
01205 n = Ndofm;
01206
01207 ni = Mp->nlman->ninr;
01208
01209 ini = Mp->nlman->niilnr;
01210
01211 zero = Mp->zero;
01212
01213 err = Mp->nlman->errnr;
01214
01215 dlambda=Mp->nlman->incrnr;
01216
01217 dlambdamin=Mp->nlman->minincrnr;
01218
01219 rb = new double [n];
01220 f = new double [n];
01221 fb = new double [n];
01222 fi = new double [n];
01223 dr = new double [n];
01224 memset (rb,0,n*sizeof(double));
01225 memset (f,0,n*sizeof(double));
01226 memset (fi,0,n*sizeof(double));
01227 memset (dr,0,n*sizeof(double));
01228
01229
01230 r = Lsrs->give_lhs (lcid);
01231 fp = Lsrs->give_rhs (lcid*2);
01232 fc = Lsrs->give_rhs (lcid*2+1);
01233 lambda=0.0;
01234
01235 i=0;
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251 for (i=0;i<ni;i++){
01252
01253
01254 for (j=0;j<n;j++){
01255 rb[j]=r[j];
01256 }
01257
01258 blambda=lambda;
01259
01260
01261
01262 fprintf (stdout,"\n lambda %e, dlambda %e",lambda,dlambda);
01263
01264
01265
01266
01267
01268 for (j=0;j<n;j++){
01269 f[j]=fc[j]+lambda*fp[j];
01270 fb[j]=dlambda*fp[j];
01271 }
01272
01273
01274 stiffness_matrix (lcid);
01275
01276
01277
01278 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
01279
01280 for (j=0;j<n;j++){
01281 r[j]+=dr[j];
01282 }
01283
01284
01285 internal_forces (lcid,fi);
01286
01287
01288 for (j=0;j<n;j++){
01289 fb[j]=pfp[j]+f[j]+dlambda*fp[j];
01290 }
01291 norfa=ss(fb,fb,n);
01292 for (j=0;j<n;j++){
01293 fb[j] -= fi[j];
01294 }
01295
01296 norf=ss(fb,fb,n);
01297 norf /= norfa;
01298
01299
01300
01301 if (norf<err){
01302 lambda+=dlambda;
01303 Mm->updateipval ();
01304
01305
01306
01307
01308
01309
01310
01311 print_step(lcid, i, lambda, f);
01312 print_flush();
01313 continue;
01314 }
01315
01316
01317 for (j=0;j<ini;j++){
01318
01319
01320
01321 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
01322
01323 for (k=0;k<n;k++){
01324 r[k]+=dr[k];
01325 }
01326
01327
01328 internal_forces (lcid,fi);
01329
01330
01331 for (k=0;k<n;k++){
01332 fb[k]=pfp[k]+f[k]+dlambda*fp[k]-fi[k];
01333 }
01334
01335
01336 norf=ss(fb,fb,n);
01337 norf /= norfa;
01338
01339 fprintf (stdout,"\n j=%ld norf=%e",j,norf);
01340
01341
01342 if (norf<err) break;
01343 }
01344
01345 if (j==ini || norf>err){
01346 dlambda/=2.0;
01347 if (dlambda<dlambdamin) dlambda=dlambdamin;
01348
01349 if (Mespr==1){
01350 fprintf (stdout,"\n increment must be modified (dlambda decrease) dlambda=%e norf=%e",dlambda,norf);
01351
01352 }
01353
01354 for (j=0;j<n;j++){
01355 r[j]=rb[j];
01356 }
01357 lambda=blambda;
01358
01359 }
01360 else{
01361 lambda+=dlambda;
01362 Mm->updateipval ();
01363 print_step(lcid, i, lambda, f);
01364 print_flush();
01365
01366
01367
01368
01369
01370
01371
01372 }
01373
01374 }
01375
01376 delete [] dr; delete [] fi; delete [] fb; delete [] f;
01377
01378
01379 return lambda;
01380 }