00001 #include <math.h>
00002 #include <string.h>
00003 #include "cpcsolver.h"
00004 #include "globalc.h"
00005 #include "global.h"
00006 #include "globalt.h"
00007 #include "globmat.h"
00008 #include "globmatt.h"
00009 #include "globmatc.h"
00010 #include "elemswitcht.h"
00011 #include "gmatrix.h"
00012 #include "gtopology.h"
00013 #include "mechprint.h"
00014 #include "transprint.h"
00015 #include "backupsol.h"
00016 #include "backupsolt.h"
00017 #include "npsolvert.h"
00018
00019
00020
00021
00022
00023
00024 void solve_gpcouplprob ()
00025 {
00026 long lcid;
00027
00028
00029
00030 lcid=0;
00031
00032 switch (Cp->tnlinsol){
00033 case newtonc:{
00034 newton_raphson_gparcoupl (lcid);
00035 break;
00036 }
00037 default:{
00038 print_err("unknown solver of nonlinear equation system is required",__FILE__,__LINE__,__func__);
00039 }
00040 }
00041 }
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 void newton_raphson_gparcoupl (long lcid)
00054 {
00055 switch (Tp->tprob){
00056 case growing_np_problem:{
00057
00058
00059 newton_raphson_gparcoupl_lin (lcid);
00060 break;
00061 }
00062 case growing_np_problem_nonlin:{
00063
00064
00065 newton_raphson_gparcoupl_nonlin (lcid);
00066 break;
00067 }
00068 default:{
00069 print_err("unknown TRFEL problem type is required",__FILE__,__LINE__,__func__);
00070 }
00071 }
00072 }
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 void newton_raphson_gparcoupl_lin (long lcid)
00085 {
00086 int mefel_convergence;
00087 long i,ii,j,k,nm,nt,ini,nsts,mncd,tncd,mnce,tnce;
00088 double zero,dt,dtmin,dtdef,end_time,alpha,s,prev_timem,prev_timet;
00089 double *dr,*fb,*r,*d,*f,*fp,*fi,*fl,*p,*lhst,*tdlhst,*rhst,*lhsb,*lhstb,*tdlhstb;
00090 double norf, norfa, err;
00091 matrix lsm_a(3,3);
00092 vector lsm_r(3), lsm_l(3);
00093
00094
00095 ini = Mp->niilnr;
00096
00097 err = Mp->errnr;
00098
00099 nm=Ndofm;
00100
00101 nt=Ndoft;
00102
00103
00104
00105
00106
00107
00108 lhst=Lsrst->give_lhs (lcid);
00109
00110 tdlhst=Lsrst->give_tdlhs (lcid);
00111
00112 rhst=Lsrst->give_rhs (lcid);
00113
00114
00115 d = new double [nt];
00116 nullv (d,nt);
00117
00118 p = new double [nt];
00119 nullv (p,nt);
00120
00121 lhstb = new double [nt];
00122 nullv (lhstb,nt);
00123
00124 tdlhstb = new double [nt];
00125 nullv (tdlhstb,nt);
00126
00127
00128
00129
00130
00131
00132 r = Lsrs->give_lhs (lcid);
00133
00134 f = Lsrs->give_rhs (lcid);
00135
00136
00137 dr = new double [nm];
00138 nullv (dr,nm);
00139
00140 fp = new double [nm];
00141 nullv (fp,nm);
00142
00143 fl = new double [nm];
00144 nullv (fl,nm);
00145
00146 fi = new double [nm];
00147 nullv (fi,nm);
00148
00149 fb = new double [nm];
00150 nullv (fb,nm);
00151
00152 lhsb = new double [nm];
00153 nullv (lhsb,nm);
00154
00155
00156 alpha=Tp->alpha;
00157 zero=Tp->zero;
00158
00159
00160
00161 Cp->time=Cp->timecon.starttime ();
00162
00163 dt=Cp->timecon.initialtimeincr ();
00164
00165 end_time = Cp->timecon.endtime ();
00166
00167
00168 Mp->time = Cp->time;
00169 Mp->timecon.take_values (Cp->timecon);
00170
00171 dtmin=Mp->timecon.dtmin;
00172
00173 prev_timem = Mp->time;
00174
00175 Tp->time = Cp->time;
00176 Tp->timecont.take_values (Cp->timecon);
00177
00178 prev_timet = Tp->time;
00179
00180
00181 approximation_inittemper();
00182
00183
00184 approximation();
00185
00186 actual_previous_nodval();
00187
00188
00189
00190
00191 i=0;
00192 ii=0;
00193
00194 nsts=0;
00195
00196
00197 if (Mp->hdbcont.restore_stat()){
00198 if (Mespr==1)
00199 fprintf (stdout,"\n Reading of backup file for MEFEL\n");
00200 solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm);
00201
00202
00203 Mm->initmaterialmodels(lcid);
00204 print_init(-1, "at");
00205
00206
00207 }
00208 else
00209 {
00210
00211 Mm->initmaterialmodels(lcid);
00212
00213 internal_forces(lcid, fp);
00214 print_init(-1, "wt");
00215 print_step(lcid, i, Mp->time, f);
00216 print_flush();
00217 }
00218 if (Tp->hdbcont.restore_stat()){
00219 if (Mesprt==1)
00220 fprintf (stdout,"\n Reading of backup file for TRFEL\n");
00221 solvert_restore (lhst,tdlhst,f,i,Cp->time,dt,Cp->timecon,Ndoft);
00222
00223 Tm->initmaterialmodels();
00224 compute_req_valt (0);
00225 print_initt(-1, "at");
00226
00227
00228 }
00229 else{
00230
00231 Tm->initmaterialmodels();
00232 compute_req_valt (0);
00233 print_initt(-1, "wt");
00234 print_stept(0,i,Tp->time,NULL);
00235 print_flusht();
00236 }
00237
00238
00239 mefel_convergence = 1;
00240
00241
00242 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00243 Mt->save_nodval(lcid);
00244
00245
00246
00247
00248
00249 do{
00250 if(mefel_convergence == 0){
00251 Cp->timecon.oldtime ();
00252 Tp->timecont.oldtime ();
00253 if (Mesprc==1) fprintf (stdout,"\n\n --------------------------------------------------------------------------------------------");
00254 if (Mesprc==1) fprintf (stdout,"\n REDUCED TIME STEP, METR TIME STEP = %ld : METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
00255 if (Mesprc==1) fprintf (stdout,"\n --------------------------------------------------------------------------------------------");
00256 i--;
00257 }
00258
00259
00260
00261 Cp->time = Cp->timecon.newtime (dt);
00262 Tp->time = Tp->timecont.newtime (dt);
00263 Tp->time = Cp->time;
00264 Tp->timecont.take_values (Cp->timecon);
00265 Tp->timecont = Cp->timecon;
00266 i++;
00267
00268 if (Cp->time >= 259200.0)
00269 fprintf (stdout,"\n trfel tady !!!");
00270
00271 if(mefel_convergence == 1){
00272 if (Mesprc==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------------------");
00273 if (Mesprc==1) fprintf (stdout,"\n NEW METR TIME STEP = %ld: METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
00274 if (Mesprc==1) fprintf (stdout,"\n -------------------------------------------------------------------------------------------");
00275 }
00276 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
00277 if (Mesprt==1) fprintf (stdout,"\n TRFEL TIME STEP = %ld, TRFEL TIME = %e, trfel time increment = %e",i,Tp->time,dt);
00278 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
00279
00280
00281
00282
00283 tnce = Gtt->search_newelem (Tp->time,prev_timet);
00284
00285
00286
00287
00288 if (tnce>0){
00289 Tt->initial_nodval ();
00290 }
00291
00292
00293
00294
00295 tncd = Gtt->search_newdofs (Tp->time,prev_timet);
00296
00297
00298 Gtt->update_elem (Tp->time);
00299
00300
00301 Gtt->update_nodes ();
00302
00303
00304 Gtt->update_dofs (Tp->time);
00305
00306
00307 Ndoft = Gtt->codenum_generation (Outt);
00308 nt=Ndoft;
00309
00310
00311 if (tncd>0){
00312 if (Kmat != NULL){
00313 delete Kmat;
00314 Kmat=NULL;
00315 }
00316
00317 if (Cmat != NULL){
00318 delete Cmat;
00319 Cmat=NULL;
00320 }
00321 }
00322
00323 if (tnce!=0 || tncd!=0){
00324
00325 Tt->lhs_restore (lhst,Lsrst->lhsi,tdlhst);
00326 }
00327
00328 Tm->updateipval ();
00329
00330
00331 conductivity_matrix (lcid);
00332
00333
00334 capacity_matrix (lcid);
00335
00336
00337 for (j=0;j<nt;j++){
00338 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
00339 }
00340
00341
00342 trfel_right_hand_side (0,rhst,nt);
00343
00344
00345 Kmat->gmxv (p,d);
00346
00347
00348
00349 Kmat->scalgm (dt*alpha);
00350 Kmat->addgm (1.0,*Cmat);
00351
00352 for (j=0;j<nt;j++){
00353 rhst[j] = rhst[j] - d[j];
00354 d[j]=tdlhst[j];
00355 }
00356
00357
00358 Tp->ssle->solve_system (Gtt,Kmat,tdlhst,rhst,Outt);
00359
00360
00361 for (j=0;j<nt;j++){
00362 lhstb[j]=lhst[j];
00363 tdlhstb[j]=tdlhst[j];
00364 }
00365
00366
00367 for (j=0;j<nt;j++){
00368 s=(1.0-alpha)*d[j]+alpha*tdlhst[j];
00369 lhst[j]+=dt*s;
00370 }
00371
00372
00373 solution_correction ();
00374
00375 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00376
00377
00378 approximation ();
00379
00380 actual_previous_nodval ();
00381
00382 if (Cp->fcsolv == fullnewtonc){
00383
00384
00385 if (Kmat != NULL){
00386 delete Kmat;
00387 Kmat=NULL;
00388 }
00389
00390 if (Cmat != NULL){
00391 delete Cmat;
00392 Cmat=NULL;
00393 }
00394 }
00395
00396
00397
00398
00399 mefel_convergence = 1;
00400
00401 Mp->time = Mp->timecon.newtime (dt);
00402
00403 if (Mp->time <= Tp->time){
00404
00405 if (Cp->time >= 259200.0)
00406 fprintf (stdout,"\n mefel tady !!!");
00407
00408
00409 ii++;
00410
00411
00412 Mp->strainstate=0;
00413
00414
00415 Mp->stressstate=0;
00416
00417
00418 Mp->otherstate=0;
00419
00420
00421 for (j=0;j<nm;j++){
00422 lhsb[j]=r[j];
00423 }
00424
00425 if (Mespr==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------");
00426 if (Mespr==1) fprintf (stdout,"\n MEFEL TIME STEP = %ld, MEFEL TIME = %e, mefel time increment = %e",ii,Mp->time,dt);
00427 if (Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------------------------");
00428
00429
00430
00431
00432 mnce = Gtm->search_newelem (Mp->time,prev_timem);
00433
00434 if (mnce!=0){
00435
00436
00437 Mt->initial_displ(lcid);
00438 }
00439
00440
00441
00442
00443 mncd = Gtm->search_newdofs (Mp->time,prev_timem);
00444
00445
00446 Gtm->update_elem (Mp->time);
00447
00448
00449 Gtm->update_nodes ();
00450
00451
00452 Gtm->update_dofs (Mp->time);
00453
00454
00455 Ndofm = Gtm->codenum_generation (Out);
00456 nm=Ndofm;
00457
00458
00459 if (mncd>0){
00460 if (Smat != NULL){
00461 delete Smat;
00462 Smat=NULL;
00463 }
00464 }
00465
00466
00467 if (mnce!=0 || mncd!=0){
00468 Mt->restore_nodval (r,nm);
00469
00470 }
00471
00472
00473 approximation_temper ();
00474
00475 approximation_humid ();
00476
00477
00478 mefel_right_hand_side (lcid,f,fl);
00479
00480
00481
00482 Mb->comp_sum (f);
00483 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fx=%le, fy=%le, fz=%le",
00484 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00485 fflush(Out);
00486
00487
00488 incr_internal_forces (lcid,fi);
00489
00490
00491 for (j=0;j<nm;j++){
00492 fb[j]=f[j]-fp[j]+fi[j];
00493 }
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505 if ((Mp->nlman->stmat==tangent_stiff) || (ii == 1) || (Smat==NULL))
00506 stiffness_matrix (lcid);
00507
00508
00509 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00510
00511
00512 for (j=0;j<nm;j++){
00513 r[j]+=dr[j];
00514 }
00515
00516
00517 internal_forces (lcid,fi);
00518
00519
00520 norfa=ss(f,f,nm);
00521
00522 for (j=0;j<nm;j++){
00523 fb[j] = fl[j] - fi[j];
00524 }
00525
00526 norf=ss(fb,fb,nm);
00527 if (norfa > 0.0)
00528 norf /= norfa;
00529
00530 if (Mespr==1) fprintf (stdout,"\n\n Norf before inner iteration loop, norf = %e\n\n\n",norf);
00531 j = 0;
00532 if (norf<err){
00533 nsts++;
00534 }
00535 else
00536 {
00537
00538
00539 fillm(0.0, lsm_a);
00540 fillv(0.0, lsm_r);
00541 for (j=0;j<ini;j++)
00542 {
00543 if((Cp->fcsolv == fullnewtonc) ||
00544 ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0)))
00545 {
00546
00547 stiffness_matrix (lcid);
00548 }
00549
00550 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00551
00552
00553 for (k=0;k<nm;k++)
00554 r[k]+=dr[k];
00555
00556
00557 internal_forces (lcid,fi);
00558
00559
00560 for (k=0;k<nm;k++)
00561 fb[k]= fl[k] - fi[k];
00562
00563
00564 norf=ss(fb,fb,nm);
00565 if (norfa > 0.0)
00566 norf /= norfa;
00567
00568
00569 if (Mespr==1)
00570 fprintf (stdout,"\n Inner loop number j = %ld, norf = %e \n",j,norf);
00571
00572 if (norf<err) break;
00573
00574
00575 if (j > 1)
00576 {
00577
00578 if (lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero, 1) > 0.0)
00579 {
00580 if (Mespr==1)
00581 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00582 break;
00583 }
00584 }
00585 else
00586 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero, 0);
00587 }
00588 }
00589
00590 if(Cp->fcsolv == fullnewtonc){
00591
00592 if (Smat != NULL){
00593 delete Smat;
00594 Smat=NULL;
00595 }
00596 }
00597
00598 if (j==ini || norf>err)
00599 {
00600
00601 if (Mespr==1) fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00602
00603 mefel_convergence = 0;
00604
00605
00606 for (j=0;j<nm;j++)
00607 r[j]=lhsb[j];
00608
00609
00610
00611 dt/=2.0;
00612 Mp->timecon.oldtime ();
00613 ii--;
00614
00615 if (Mespr==1) fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00616 if (dt<dtmin){
00617 if (Mespr==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00618 if (Mespr==1) fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
00619 compute_req_val (lcid);
00620 print_step_forced(lcid, i, Mp->time, f);
00621 print_flush();
00622 break;
00623 }
00624 }
00625 else
00626 {
00627
00628 if (Mespr==1) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00629 mefel_convergence = 1;
00630
00631
00632 for (j=0;j<nm;j++)
00633 fp[j]=f[j];
00634
00635
00636
00637 Mt->save_nodval(lcid);
00638
00639 prev_timem = Mp->time;
00640
00641 dtdef = Mp->timecon.actualforwtimeincr ();
00642 if (nsts==2){
00643 dt*=2.0;
00644 nsts=0;
00645 if (Mespr==1) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00646 }
00647 if (dt>dtdef)
00648 dt=dtdef;
00649
00650
00651 Mb->comp_sum (fi);
00652 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fix=%le, fiy=%le, fiz=%le",
00653 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00654 fflush(Out);
00655
00656
00657 Mm->updateipval();
00658 compute_req_val (lcid);
00659 print_step(lcid, ii, Mp->time, f);
00660 print_flush();
00661
00662 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.hdbtype))
00663 {
00664 if (Mespr==1)
00665 fprintf (stdout,"\n Creating backup file for MEFEL\n");
00666 solver_save (r,f,ii,Mp->time,dt,&Mp->timecon,Ndofm);
00667 }
00668 }
00669 }
00670
00671 if(mefel_convergence == 0){
00672
00673 for (j=0;j<nt;j++){
00674 lhst[j]=lhstb[j];
00675 tdlhst[j]=tdlhstb[j];
00676 }
00677 }
00678 else{
00679 print_stept(0,i,Tp->time,NULL);
00680 print_flusht();
00681 if ((Tp->timecont.isitimptime()==1) && Tp->hdbcont.save_stat())
00682 {
00683 if (Mesprt==1)
00684 fprintf (stdout,"\n Creating backup file for TRFEL\n");
00685 solvert_save (lhst,tdlhst,rhst,i,Cp->time,dt,Cp->timecon,Ndoft);
00686 }
00687 }
00688 } while(Cp->time<=end_time);
00689
00690 print_close ();
00691 print_closet ();
00692
00693 delete [] fi;
00694 delete [] fb;
00695 delete [] fp;
00696 delete [] dr;
00697 delete [] fl;
00698 delete [] lhsb;
00699
00700 delete [] p;
00701 delete [] d;
00702
00703 delete [] lhstb;
00704 delete [] tdlhstb;
00705 }
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719 void newton_raphson_gparcoupl_nonlin (long lcid)
00720 {
00721 int mefel_convergence;
00722 long i,ii,j,k,nm,nt,ini,init,nsts,mncd,tncd,mnce,mnae,tnce;
00723 long *ifnm, nifnm;
00724 double zero,dt,dtmin,dtdef,end_time,alpha,prev_timem,prev_timet;
00725 double *dr,*fb,*r,*d,*f,*fp,*fi,*fl,*p,*lhst,*tdlhst,*rhst,*lhsb,*lhstb,*tdlhstb;
00726 double *fbt,*fit,*lhst_last;
00727
00728 double norf_last;
00729 double norf, norfa,err,errt;
00730 matrix lsm_a(3,3);
00731 vector lsm_r(3), lsm_l(3);
00732
00733 ifnm = new long[Mt->nn];
00734 memset(ifnm, 0, sizeof(*ifnm)*Mt->nn);
00735 Mp->nlman->divc_step = new double[10];
00736 Mp->nlman->divc_err = new double[10];
00737
00738
00739 ini = Mp->niilnr;
00740
00741 err = Mp->errnr;
00742
00743 init = Tp->nii;
00744
00745 errt = Tp->err;
00746
00747 nm=Ndofm;
00748
00749 nt=Ndoft;
00750
00751
00752
00753
00754
00755 lhst=Lsrst->give_lhs (lcid);
00756
00757 tdlhst=Lsrst->give_tdlhs (lcid);
00758
00759 rhst=Lsrst->give_rhs (lcid);
00760
00761
00762 d = new double [nt];
00763 nullv (d,nt);
00764
00765 p = new double [nt];
00766 nullv (p,nt);
00767 fbt = new double [nt];
00768 nullv (fbt,nt);
00769 fit = new double [nt];
00770 nullv (fit,nt);
00771 lhst_last = new double [nt];
00772
00773
00774 lhstb = new double [nt];
00775 nullv (lhstb,nt);
00776
00777 tdlhstb = new double [nt];
00778 nullv (tdlhstb,nt);
00779
00780
00781
00782
00783
00784
00785 r = Lsrs->give_lhs (lcid);
00786
00787 f = Lsrs->give_rhs (lcid);
00788
00789
00790 dr = new double [nm];
00791 nullv (dr,nm);
00792
00793 fp = new double [nm];
00794 nullv (fp,nm);
00795
00796 fl = new double [nm];
00797 nullv (fl,nm);
00798
00799 fi = new double [nm];
00800 nullv (fi,nm);
00801
00802 fb = new double [nm];
00803 nullv (fb,nm);
00804
00805 lhsb = new double [nm];
00806 nullv (lhsb,nm);
00807
00808
00809 alpha=Tp->alpha;
00810 zero=Tp->zero;
00811
00812
00813
00814 Cp->time=Cp->timecon.starttime ();
00815
00816 dt=Cp->timecon.initialtimeincr ();
00817
00818 end_time = Cp->timecon.endtime ();
00819
00820
00821 Mp->time = Cp->time;
00822 Mp->timecon.take_values (Cp->timecon);
00823
00824 dtmin=Mp->timecon.dtmin;
00825
00826 prev_timem = Mp->time;
00827
00828 Tp->time = Cp->time;
00829 Tp->timecont.take_values (Cp->timecon);
00830
00831 prev_timet = Tp->time;
00832
00833
00834 approximation_inittemper ();
00835
00836
00837 approximation ();
00838
00839 actual_previous_nodval ();
00840
00841
00842
00843
00844 i=0;
00845 ii=0;
00846
00847 nsts=0;
00848
00849
00850 if (Mp->hdbcont.restore_stat()){
00851 if (Mespr==1)
00852 fprintf (stdout,"\n Reading of backup file for MEFEL\n");
00853 solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm);
00854
00855 prev_timem = Mp->time;
00856
00857
00858 copyv(r, lhsb, nm);
00859
00860 Gtm->update_elem (Mp->time);
00861
00862 Gtm->update_auxinf();
00863
00864 Gtm->update_nodes ();
00865
00866 Gtm->update_dofs (Mp->time);
00867
00868 nm = Ndofm = Gtm->codenum_generation(Out);
00869
00870 Mt->save_nodval(lcid);
00871
00872 Mt->comreac();
00873
00874 Mt->elemprescdisp();
00875
00876 Mm->initmaterialmodels(lcid);
00877
00878 print_init(-1, "at");
00879
00880
00881 }
00882 else
00883 {
00884
00885 Mm->initmaterialmodels(lcid);
00886
00887 internal_forces(lcid, fp);
00888 print_init(-1, "wt");
00889 print_step(lcid, i, Mp->time, f);
00890 print_flush();
00891 }
00892 if (Tp->hdbcont.restore_stat()){
00893 if (Mesprt==1)
00894 fprintf (stdout,"\n Reading of backup file for TRFEL\n");
00895 solvert_restore (lhst,tdlhst,f,i,Cp->time,dt,Cp->timecon,Ndoft);
00896
00897 Tm->initmaterialmodels();
00898 compute_req_valt (0);
00899 print_initt(-1, "at");
00900
00901
00902 }
00903 else{
00904
00905 Tm->initmaterialmodels();
00906 compute_req_valt (0);
00907 print_initt(-1, "wt");
00908 print_stept(0,i,Tp->time,NULL);
00909 print_flusht();
00910 }
00911
00912
00913 mefel_convergence = 1;
00914
00915
00916 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
00917 Mt->save_nodval(lcid);
00918
00919
00920
00921
00922
00923 do{
00924 if(mefel_convergence == 0){
00925 Cp->timecon.oldtime ();
00926 Tp->timecont.oldtime ();
00927 if (Mesprc==1) fprintf (stdout,"\n\n --------------------------------------------------------------------------------------------");
00928 if (Mesprc==1) fprintf (stdout,"\n REDUCED TIME STEP, METR TIME STEP = %ld : METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
00929 if (Mesprc==1) fprintf (stdout,"\n --------------------------------------------------------------------------------------------");
00930 i--;
00931 }
00932
00933
00934
00935 Cp->time = Cp->timecon.newtime (dt);
00936 Tp->time = Tp->timecont.newtime (dt);
00937 Tp->time = Cp->time;
00938 Tp->timecont.take_values (Cp->timecon);
00939 i++;
00940
00941 if(mefel_convergence == 1){
00942 if (Mesprc==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------------------");
00943 if (Mesprc==1) fprintf (stdout,"\n NEW METR TIME STEP = %ld: METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
00944 if (Mesprc==1) fprintf (stdout,"\n -------------------------------------------------------------------------------------------");
00945 }
00946 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
00947 if (Mesprt==1) fprintf (stdout,"\n TRFEL TIME STEP = %ld, TRFEL TIME = %e, trfel time increment = %e",i,Tp->time,dt);
00948 if (Mesprt==1) fprintf (stdout,"\n ------------------------------------------------------------------------------------");
00949
00950
00951
00952
00953 tnce = Gtt->search_newelem (Tp->time,prev_timet);
00954
00955
00956
00957
00958 if (tnce>0){
00959 Tt->initial_nodval ();
00960 }
00961
00962
00963
00964
00965 tncd = Gtt->search_newdofs (Tp->time,prev_timet);
00966
00967
00968 Gtt->update_elem (Tp->time);
00969
00970
00971 Gtt->update_nodes ();
00972
00973
00974 Gtt->update_dofs (Tp->time);
00975
00976
00977
00978 nt=Ndoft;
00979
00980
00981 if (tncd>0){
00982 if (Kmat != NULL){
00983 delete Kmat;
00984 Kmat=NULL;
00985 }
00986
00987 if (Cmat != NULL){
00988 delete Cmat;
00989 Cmat=NULL;
00990 }
00991 }
00992
00993 if (tnce!=0 || tncd!=0){
00994
00995 Tt->lhs_restore (lhst,Lsrst->lhsi,tdlhst);
00996 }
00997
00998 Tm->updateipval ();
00999
01000
01001 conductivity_matrix (lcid);
01002
01003
01004 capacity_matrix (lcid);
01005
01006
01007 for (j=0;j<nt;j++){
01008 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
01009 }
01010
01011
01012 trfel_right_hand_side (0,rhst,nt);
01013
01014
01015 Kmat->gmxv (p,d);
01016
01017
01018
01019 Kmat->scalgm (dt*alpha);
01020 Kmat->addgm (1.0,*Cmat);
01021
01022 for (j=0;j<nt;j++){
01023
01024 rhst[j] = rhst[j] - d[j];
01025
01026
01027 }
01028
01029
01030 Tp->ssle->solve_system (Gtt,Kmat,tdlhst,rhst,Outt);
01031
01032
01033 for (j=0;j<nt;j++){
01034 lhstb[j]=lhst[j];
01035 tdlhstb[j]=tdlhst[j];
01036 }
01037
01038
01039 for (j=0;j<nt;j++){
01040
01041
01042
01043
01044 lhst[j] = p[j] + alpha*dt*tdlhst[j];
01045 }
01046
01047
01048 solution_correction ();
01049
01050 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01051
01052
01053 approximation ();
01054
01055 norf_last = 1.0e20;
01056
01057 for (j=0;j<init;j++){
01058
01059 solution_correction ();
01060
01061 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01062
01063
01064 approximation ();
01065
01066
01067 if (Tp->trsolv == fullnewtont){
01068
01069
01070 capacity_matrix (0);
01071
01072
01073 conductivity_matrix (0);
01074
01075
01076 Kmat->gmxv (p,d);
01077
01078
01079
01080 Kmat->scalgm (dt*alpha);
01081 Kmat->addgm (1.0,*Cmat);
01082 }
01083
01084 if (Tp->trestype==lrhst){
01085
01086 trfel_right_hand_side (0,rhst,nt);
01087
01088 Kmat->gmxv (tdlhst,fit);
01089
01090 for (k=0;k<nt;k++){
01091 fbt[k] = rhst[k] - d[k] - fit[k];
01092 }
01093 }
01094
01095 if (Tp->trestype==fluxest){
01096
01097 internal_fluxes (fit,nt);
01098
01099 for (k=0;k<nt;k++){
01100 fbt[k]=fit[k];
01101 }
01102 }
01103
01104 norf = ss (fbt,fbt,nt);
01105
01106 if (Mesprt==1) fprintf (stdout,"\n TRFEL inner loop number j = %ld, norf = %e",j,norf);
01107
01108 if (norf<errt) break;
01109
01110 Tp->ssle->solve_system (Gtt,Kmat,d,fbt,Outt);
01111
01112 for (k=0;k<nt;k++){
01113 tdlhst[k]+=d[k];
01114 lhst[k]+=alpha*dt*d[k];
01115 }
01116
01117
01118 solution_correction ();
01119
01120 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01121
01122
01123 approximation ();
01124
01125
01126
01127 if (Tp->convergcontrolt==yes){
01128 if (norf > norf_last){
01129 for (k=0;k<nt;k++){
01130 lhst[k] = lhst_last[k];
01131 }
01132
01133
01134 solution_correction ();
01135
01136 Tt->lhs_save (lhst,Lsrst->lhsi,tdlhst);
01137
01138
01139 approximation ();
01140
01141 if (Mesprt==1) fprintf (stdout,"\n\nTRFEL convergence control: inner iteration skipped %ld error %e\n\n",j,norf);
01142
01143 break;
01144 }
01145 for (k=0;k<nt;k++){
01146 lhst_last[k] = lhst[k];
01147 }
01148 norf_last = norf;
01149 }
01150 }
01151
01152 if (Cp->fcsolv == fullnewtonc){
01153
01154
01155 if (Kmat != NULL){
01156 delete Kmat;
01157 Kmat=NULL;
01158 }
01159
01160 if (Cmat != NULL){
01161 delete Cmat;
01162 Cmat=NULL;
01163 }
01164 }
01165
01166
01167
01168
01169 mefel_convergence = 1;
01170
01171 Mp->time = Mp->timecon.newtime (dt);
01172 if (prev_timem == Mp->timecon.starttime())
01173 prev_timem = Mp->time;
01174
01175 if (Mp->time <= Tp->time){
01176
01177 ii++;
01178 Mp->istep = i;
01179
01180
01181 Mp->strainstate=0;
01182
01183
01184 Mp->stressstate=0;
01185
01186
01187 Mp->otherstate=0;
01188
01189 if (Mespr==1) fprintf (stdout,"\n\n -----------------------------------------------------------------------------------");
01190 if (Mespr==1) fprintf (stdout,"\n MEFEL TIME STEP = %ld, MEFEL TIME = %e, mefel time increment = %e",ii,Mp->time,dt);
01191 if (Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------------------------");
01192
01193
01194
01195 mnce = Gtm->search_changed_elem(Mp->time, mnae);
01196
01197
01198 mncd = Gtm->search_changed_dofs(Mp->time,prev_timem);
01199
01200 nifnm = 0;
01201 if (mnce!=0){
01202
01203
01204 nifnm = Gtm->search_iface_nodes(ifnm);
01205 }
01206
01207 if ((mnce!=0) || (mncd!=0))
01208 {
01209
01210
01211 Mt->save_node_inidispl(lcid, Mp->time, prev_timem);
01212
01213
01214 Gtm->switch_new_elem();
01215
01216
01217 Gtm->update_nodes();
01218
01219
01220
01221 Gtm->update_active_dofs(Mp->time);
01222
01223
01224 Ndofm = Gtm->codenum_generation(Out);
01225
01226
01227
01228 Mt->restore_nodval(r, Ndofm);
01229
01230
01231 Mm->initmaterialmodels(lcid);
01232
01233
01234
01235
01236
01237
01238 if (nifnm)
01239 {
01240 if (Mp->comp_inidispl)
01241 {
01242
01243
01244 stress_initdispl(lcid);
01245
01246
01247
01248 Gtm->update_active_dofs(Mp->time);
01249
01250 Gtm->clear_intf_dofs(ifnm);
01251
01252
01253 Ndofm = Gtm->codenum_generation(Out);
01254
01255
01256 Mp->strcomp = 0;
01257 internal_forces(lcid,fb);
01258 cmulv(-1.0, fb, Ndofm);
01259 Mp->strcomp = 1;
01260
01261
01262
01263 delete Smat;
01264 Smat=NULL;
01265 stiffness_matrix (lcid);
01266 nullv(r, Ndofm);
01267
01268
01269
01270 Mp->ssle->solve_system (Gtm,Smat,r,fb,Out);
01271
01272
01273
01274
01275
01276
01277 Mt->save_elem_inidispl(lcid, ifnm);
01278
01279
01280
01281 Gtm->remove_nodes(ifnm);
01282
01283
01284
01285 Mt->save_nodval(lcid);
01286 }
01287 else
01288 {
01289
01290
01291
01292
01293
01294
01295 Mt->save_elem_inidispl(lcid, ifnm);
01296 }
01297 }
01298
01299 Gtm->update_elem(Mp->time);
01300
01301
01302 Gtm->update_nodes();
01303
01304
01305 Gtm->update_dofs (Mp->time);
01306
01307
01308 nm = Ndofm = Gtm->codenum_generation(Out);
01309
01310
01311 delete Smat;
01312 Smat=NULL;
01313
01314
01315 Mt->restore_nodval(r,nm);
01316
01317
01318
01319
01320 if (i > 1)
01321 Mt->restore_nodforce(fp, nm);
01322
01323
01324 Mt->comreac();
01325
01326 Mt->elemprescdisp();
01327
01328 Mt->clean_ip_new_elem();
01329 }
01330
01331
01332 approximation_temper ();
01333
01334 approximation_humid ();
01335
01336
01337 mefel_right_hand_side (lcid,f,fl);
01338
01339
01340
01341 Mb->comp_sum (f);
01342 Mb->comp_sum_react();
01343 Mb->comp_sum_pdreact();
01344 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fx=%le, fy=%le, fz=%le",
01345 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
01346 fflush(Out);
01347
01348 if ((mnce - mnae) > 0)
01349 {
01350
01351
01352
01353
01354 Mp->strcomp = 0;
01355 internal_forces(lcid,fi);
01356 Mp->strcomp = 1;
01357
01358
01359 subv(fp, fi, fb, nm);
01360 }
01361 else
01362 nullv(fb, nm);
01363
01364
01365 incr_internal_forces (lcid,fi);
01366
01367
01368 for (j=0;j<nm;j++){
01369 fb[j]+=f[j]-fp[j]+fi[j];
01370 }
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382 if ((Mp->nlman->stmat==tangent_stiff) || (ii == 1) || (Smat==NULL))
01383 stiffness_matrix (lcid);
01384
01385
01386 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
01387
01388
01389 addv(r, dr, nm);
01390
01391
01392 internal_forces (lcid,fi);
01393
01394
01395 norfa=ss(f,f,nm);
01396
01397
01398 subv(fl, fi, fb, nm);
01399
01400
01401 norf=ss(fb,fb,nm);
01402 if (norfa > 0.0)
01403 norf /= norfa;
01404
01405 if (Mespr==1) fprintf (stdout,"\n\n Norf before inner iteration loop norf = %e\n\n\n",norf);
01406 j = 0;
01407 if (norf<err){
01408
01409 nsts++;
01410 }
01411 else
01412 {
01413
01414
01415 fillm(0.0, lsm_a);
01416 fillv(0.0, lsm_r);
01417 for (j=0;j<ini;j++)
01418 {
01419 if((Cp->fcsolv == fullnewtonc) ||
01420 ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0)))
01421 {
01422
01423 stiffness_matrix (lcid);
01424 }
01425
01426 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
01427
01428
01429 for (k=0;k<nm;k++)
01430 r[k]+=dr[k];
01431
01432
01433 internal_forces (lcid,fi);
01434
01435
01436 subv(fl, fi, fb, nm);
01437
01438
01439 norf=ss(fb,fb,nm);
01440 if (norfa > 0.0)
01441 norf /= norfa;
01442
01443
01444 if (Mespr==1)
01445 fprintf (stdout,"\n Inner loop j = %ld norf = %e\n",j,norf);
01446
01447 if (norf<err) break;
01448
01449
01450 if (j > 1)
01451 {
01452 if (lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero, 1) > 0.0)
01453 {
01454 if (Mespr==1)
01455 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
01456 break;
01457 }
01458 }
01459 else
01460 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero, 0);
01461 }
01462 }
01463
01464 if(Cp->fcsolv == fullnewtonc){
01465
01466 if (Smat != NULL){
01467 delete Smat;
01468 Smat=NULL;
01469 }
01470 }
01471
01472 if (j==ini || norf>err)
01473 {
01474
01475 if (Mespr==1) fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
01476
01477 mefel_convergence = 0;
01478
01479
01480
01481 dt/=2.0;
01482 Mp->timecon.oldtime ();
01483 ii--;
01484
01485
01486
01487
01488 if ((mnce!=0) || (mncd!=0))
01489 {
01490
01491 Gtm->update_elem(prev_timem);
01492
01493 Gtm->update_nodes();
01494
01495 Gtm->update_dofs (prev_timem);
01496
01497 nm = Ndofm = Gtm->codenum_generation(Out);
01498 }
01499
01500
01501 for (j=0;j<nm;j++)
01502 r[j]=lhsb[j];
01503
01504 if (Mespr==1) fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
01505 if (dt<dtmin){
01506 if (Mespr==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
01507 if (Mespr==1) fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
01508 compute_req_val (lcid);
01509 print_step_forced(lcid, i, Mp->time, f);
01510 print_flush();
01511 break;
01512 }
01513 }
01514 else
01515 {
01516
01517 if (Mespr==1) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
01518 mefel_convergence = 1;
01519
01520
01521
01522 copyv(r, lhsb, nm);
01523
01524
01525
01526 copyv(f, fp, nm);
01527
01528
01529
01530
01531
01532
01533
01534 Mt->save_nodval(lcid);
01535
01536 Mt->save_nodforce(f);
01537
01538
01539 Gtm->update_auxinf();
01540 prev_timem = Mp->time;
01541
01542 dtdef = Mp->timecon.actualforwtimeincr ();
01543 if (nsts==2){
01544 dt*=2.0;
01545 nsts=0;
01546 if (Mespr==1) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
01547 }
01548 if (dt>dtdef)
01549 dt=dtdef;
01550
01551
01552 Mb->comp_sum (fi);
01553 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fix=%le, fiy=%le, fiz=%le",
01554 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
01555 fflush(Out);
01556
01557
01558 Mm->updateipval();
01559 compute_req_val (lcid);
01560 print_step(lcid, ii, Mp->time, f);
01561 print_flush();
01562
01563 Mb->comp_sum_react();
01564 Mb->comp_sum_pdreact();
01565 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], rx=%le, ry=%le, rz=%le",
01566 i, Mp->time, Mp->time/86400, Mb->reactsumcomp[0], Mb->reactsumcomp[1], Mb->reactsumcomp[2]);
01567 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], pdx=%le, pdy=%le, pdz=%le",
01568 i, Mp->time, Mp->time/86400, Mb->pd_reactsumcomp[0], Mb->pd_reactsumcomp[1], Mb->pd_reactsumcomp[2]);
01569 fflush(Out);
01570
01571 if ((Mp->timecon.isitimptime ()==1) && (Mp->hdbcont.hdbtype))
01572 {
01573 if (Mespr==1)
01574 fprintf (stdout,"\n Creating backup file for MEFEL\n");
01575 solver_save (r,f,ii,Mp->time,dt,&Mp->timecon,Ndofm);
01576 }
01577 }
01578 }
01579 if(mefel_convergence == 0){
01580
01581 for (j=0;j<nt;j++){
01582 lhst[j]=lhstb[j];
01583 tdlhst[j]=tdlhstb[j];
01584 }
01585 }
01586 else{
01587 print_stept(0,i,Tp->time,NULL);
01588 print_flusht();
01589 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat()){
01590 if (Mesprt==1)
01591 fprintf (stdout,"\n Creating backup file for TRFEL\n");
01592 solvert_save (lhst,tdlhst,rhst,i,Cp->time,dt,Cp->timecon,Ndoft);
01593 }
01594 }
01595 }while(Cp->time<end_time);
01596
01597 print_close ();
01598 print_closet ();
01599
01600 delete [] fi;
01601 delete [] fb;
01602 delete [] fp;
01603 delete [] dr;
01604 delete [] fl;
01605 delete [] lhsb;
01606 delete [] ifnm;
01607
01608 delete [] p;
01609 delete [] d;
01610
01611 delete [] fit;
01612 delete [] fbt;
01613 delete [] lhst_last;
01614
01615 delete [] lhstb;
01616 delete [] tdlhstb;
01617 }