00001 #include "nnpsolvert.h"
00002 #include "npsolvert.h"
00003 #include "globalt.h"
00004 #include "globmatt.h"
00005 #include "elemswitcht.h"
00006 #include "transprint.h"
00007 #include "backupsolt.h"
00008 #include <string.h>
00009
00010 void solve_nonlinear_nonstationary_problem ()
00011 {
00012
00013
00014 nonlinear_nonstat_solv ();
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 }
00028
00029
00030
00031
00032
00033
00034
00035 void nonlinear_nonstat_solv ()
00036 {
00037 long i,j,k,nt,ini;
00038
00039 double zero,dt,end_time,alpha;
00040 double *p,*d,*lhst,*tdlhst,*rhst;
00041 double *fbt,*fit,*lhst_last;
00042
00043 double norf_last;
00044 double norf,*err,*thresh;
00045
00046
00047 ini = Tp->nii;
00048
00049 err = Tp->errarr;
00050
00051 thresh=Tp->threshrhs;
00052
00053 nt=Ndoft;
00054
00055
00056 lhst=Lsrst->give_lhs (0);
00057 tdlhst=Lsrst->give_tdlhs (0);
00058
00059 rhst=Lsrst->give_rhs (0);
00060
00061
00062 d = new double [nt];
00063
00064 p = new double [nt];
00065 fbt = new double [nt];
00066 fit = new double [nt];
00067 lhst_last = new double [nt];
00068
00069
00070 nullv (lhst,nt);
00071 nullv (tdlhst,nt);
00072 nullv (d,nt);
00073 nullv (p,nt);
00074 nullv (fbt,nt);
00075 nullv (fit,nt);
00076
00077
00078 alpha=Tp->alpha;
00079
00080 zero=Tp->zero;
00081
00082
00083 Tp->time=Tp->timecont.starttime ();
00084
00085 dt=Tp->timecont.initialtimeincr ();
00086
00087 end_time = Tp->timecont.endtime ();
00088
00089
00090 Tm->initmaterialmodels();
00091
00092 approximation ();
00093
00094
00095
00096 if (Tp->nvs==1 && Tp->pnvs==1)
00097 actual_previous_nodval ();
00098
00099
00100
00101
00102 i=0;
00103
00104
00105 if (Tp->hdbcont.restore_stat()){
00106 solvert_restore (lhst,tdlhst,rhst,i,Tp->time,dt,Tp->timecont,nt);
00107 print_initt(-1, "at");
00108
00109 }
00110 else{
00111 compute_req_valt (0);
00112 print_initt(-1, "wt");
00113 print_stept(0,i,Tp->time,rhst);
00114 }
00115
00116 do{
00117
00118 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00119
00120 Tm->updateipval ();
00121 i++;
00122
00123
00124 conductivity_matrix (0);
00125
00126
00127
00128 capacity_matrix (0);
00129
00130 for (j=0;j<nt;j++){
00131 p[j] = lhst[j] + (1.0-alpha)*dt*tdlhst[j];
00132 }
00133
00134
00135 trfel_right_hand_side (0,rhst,nt);
00136
00137 fprintf (Outt,"\n\n\n\n");
00138 for (j=0;j<nt;j++){
00139 fprintf (Outt,"\n rhst %4ld %lf",j,rhst[j]);
00140 }
00141 fprintf (Outt,"\n\n\n\n");
00142
00143
00144 Kmat->gmxv (p,d);
00145
00146
00147
00148 Kmat->scalgm (dt*alpha);
00149 Kmat->addgm (1.0,*Cmat);
00150 Kmat->copygm (*Cmat);
00151
00152 for (j=0;j<nt;j++){
00153 rhst[j] = rhst[j] - d[j];
00154 }
00155
00156
00157 Tp->ssle->solve_system (Gtt,Cmat,tdlhst,rhst,Outt);
00158
00159 for (j=0;j<nt;j++){
00160
00161
00162
00163
00164 lhst[j] = p[j] + alpha*dt*tdlhst[j];
00165 }
00166
00167
00168 solution_correction ();
00169
00170 approximation ();
00171
00172
00173
00174 if (Tp->nvs==1 && Tp->pnvs==1)
00175 actual_previous_nodval ();
00176
00177 norf_last = 1.0e20;
00178
00179 for (j=0;j<ini;j++){
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 if (Tp->trsolv == fullnewtont){
00196
00197
00198 capacity_matrix (0);
00199
00200
00201 conductivity_matrix (0);
00202
00203
00204 Kmat->gmxv (p,d);
00205
00206
00207
00208 Kmat->scalgm (dt*alpha);
00209 Kmat->addgm (1.0,*Cmat);
00210 }
00211
00212 if (Tp->trestype==lrhst){
00213
00214 trfel_right_hand_side (0,rhst,nt);
00215
00216 Kmat->gmxv (tdlhst,fit);
00217
00218 for (k=0;k<nt;k++){
00219 fbt[k] = rhst[k] - d[k] - fit[k];
00220 }
00221 }
00222
00223 if (Tp->trestype==fluxest){
00224
00225 internal_fluxes (fit,nt);
00226
00227 for (k=0;k<nt;k++){
00228 fbt[k]=fit[k];
00229 }
00230 }
00231
00232 norf = ss (fbt,fbt,nt);
00233 if (Mesprt==1) fprintf (stdout,"\n inner iteration %ld error %e",j,norf);
00234 if (norf<err[0]) break;
00235
00236
00237
00238
00239
00240
00241 Tp->ssle->solve_system (Gtt,Kmat,d,fbt,Outt);
00242
00243 for (k=0;k<nt;k++){
00244 tdlhst[k]+=d[k];
00245 lhst[k]+=alpha*dt*d[k];
00246 }
00247
00248
00249 solution_correction ();
00250
00251 approximation ();
00252
00253
00254
00255 if (Tp->nvs==1 && Tp->pnvs==1)
00256 actual_previous_nodval ();
00257
00258
00259
00260 if (Tp->convergcontrolt==yes){
00261 if (norf > norf_last){
00262 for (k=0;k<nt;k++){
00263 lhst[k] = lhst_last[k];
00264 }
00265
00266
00267 solution_correction ();
00268
00269 approximation ();
00270
00271
00272
00273 if (Tp->nvs==1 && Tp->pnvs==1)
00274 actual_previous_nodval ();
00275
00276 if (Mesprt==1) fprintf (stdout,"\n\n convergence control: inner iteration skiped %ld error %e\n\n",j,norf);
00277
00278 break;
00279 }
00280 for (k=0;k<nt;k++){
00281 lhst_last[k] = lhst[k];
00282 }
00283
00284 norf_last = norf;
00285 }
00286 }
00287
00288
00289 Tp->time = Tp->timecont.newtime ();
00290 dt = Tp->timecont.actualbacktimeincr ();
00291
00292
00293 compute_req_valt (0);
00294 print_stept(0,i,Tp->time,rhst);
00295 print_flusht();
00296
00297 if ((Tp->timecont.isitimptime ()==1) && Tp->hdbcont.save_stat())
00298 solvert_save (lhst,tdlhst,rhst,i,Tp->time,dt,Tp->timecont,nt);
00299
00300 }while(Tp->time<=end_time);
00301
00302 print_closet ();
00303
00304 delete [] p;
00305 delete [] d;
00306
00307 delete [] fit;
00308 delete [] fbt;
00309 delete [] lhst_last;
00310 }
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 void nonlinear_nonstat_solv_oldd ()
00323 {
00324 long i,j,k,n,ini,stop;
00325 double dt,end_time,alpha,norfb,err;
00326 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*arhs;
00327
00328 n=Ndoft;
00329
00330 lhs = Lsrst->give_lhs (0);
00331 tdlhs = Lsrst->give_tdlhs (0);
00332 lhsi = Lsrst->lhsi;
00333 rhs = Lsrst->give_rhs (0);
00334
00335 d = new double [n];
00336 arhs = new double [n];
00337 p = new double [n];
00338 fb = new double [n];
00339 fi = new double [n];
00340
00341
00342 nullv (lhs,n);
00343 nullv (tdlhs,n);
00344 nullv (rhs,n);
00345 nullv (arhs,n);
00346 nullv (d,n);
00347 nullv (p,n);
00348 nullv (fb,n);
00349 nullv (fi,n);
00350
00351
00352 approximation ();
00353
00354 alpha=Tp->alpha;
00355 err=Tp->err;
00356 ini=Tp->nii;
00357
00358
00359 Tp->time=Tp->timecont.starttime ();
00360
00361 end_time = Tp->timecont.endtime ();
00362
00363 dt = Tp->timecont.initialtimeincr ();
00364
00365
00366
00367
00368 i=0;
00369
00370 Tm->initmaterialmodels();
00371 print_initt(-1, "wt");
00372 print_stept(0,i,Tp->time,NULL);
00373 do{
00374
00375 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00376
00377 Tm->updateipval ();
00378
00379 fprintf (Outt,"\n\n krok %ld",i);
00380 fprintf (Outt,"\n %le %le %le %le",Tm->ip[0].av[0],Tm->ip[1].av[0],Tm->ip[4].av[0],Tm->ip[5].av[0]);
00381
00382
00383 capacity_matrix (0);
00384
00385
00386 conductivity_matrix (0);
00387
00388 for (j=0;j<n;j++){
00389 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00390 }
00391
00392
00393 trfel_right_hand_side (0,rhs,n);
00394
00395
00396 Kmat->gmxv (d,p);
00397
00398
00399
00400 Kmat->scalgm (dt*alpha);
00401 Kmat->addgm (1.0,*Cmat);
00402
00403
00404 for (j=0;j<n;j++){
00405 fb[j] = rhs[j] - p[j];
00406 }
00407
00408
00409 Tp->ssle->solve_system (Gtt,Kmat,tdlhs,fb,Outt);
00410
00411
00412 for (j=0;j<n;j++){
00413 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00414 }
00415
00416 solution_correction ();
00417 approximation ();
00418 internal_fluxes (fi,n);
00419
00420
00421 for (j=0;j<n;j++){
00422 fb[j]=fi[j];
00423
00424 }
00425
00426 norfb = ss (fb,fb,n);
00427
00428 stop = norm_computation (fb,rhs,err,i,0);
00429
00430
00431
00432
00433
00434 if (stop==Tp->ntm){
00435
00436 Tp->time = Tp->timecont.newtime ();
00437 dt = Tp->timecont.actualbacktimeincr ();
00438
00439 i++;
00440 print_stept(0,i,Tp->time,NULL);
00441 print_flusht();
00442 continue;
00443
00444 }
00445
00446 for (j=0;j<ini;j++){
00447
00448
00449 Tp->ssle->solve_system (Gtt,Kmat,p,fb,Outt);
00450
00451 for (k=0;k<n;k++){
00452 tdlhs[k]+=p[k];
00453 lhs[k]+=alpha*dt*p[k];
00454 }
00455
00456 solution_correction ();
00457 approximation ();
00458 internal_fluxes (fi,n);
00459
00460
00461 for (k=0;k<n;k++){
00462 fb[k]= fi[k];
00463 }
00464
00465
00466 norfb = ss (fb,fb,n);
00467
00468
00469
00470
00471 if (stop==Tp->ntm){
00472 break;
00473 }
00474 }
00475
00476
00477 Tp->time = Tp->timecont.newtime ();
00478 dt = Tp->timecont.actualbacktimeincr ();
00479 i++;
00480
00481 print_stept(0,i,Tp->time,NULL);
00482 print_flusht();
00483
00484 }while(Tp->time<end_time);
00485
00486 delete [] fi; delete [] fb; delete [] p;
00487 delete [] arhs; delete [] d;
00488
00489 print_closet();
00490 }
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 void nonlinear_nonstat_solv_vform ()
00503 {
00504 long i,j,k,n,ini,stop;
00505 double dt,end_time,alpha,err;
00506 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*arhs;
00507
00508
00509 n=Ndoft;
00510
00511
00512 lhs = Lsrst->give_lhs (0);
00513
00514 tdlhs = Lsrst->give_tdlhs (0);
00515
00516 lhsi = Lsrst->lhsi;
00517
00518 rhs = Lsrst->give_rhs (0);
00519
00520
00521 d = new double [n];
00522 arhs = new double [n];
00523 p = new double [n];
00524 fb = new double [n];
00525 fi = new double [n];
00526
00527
00528 nullv (lhs,n);
00529 nullv (tdlhs,n);
00530 nullv (rhs,n);
00531 nullv (arhs,n);
00532 nullv (d,n);
00533 nullv (p,n);
00534 nullv (fb,n);
00535 nullv (fi,n);
00536
00537
00538 approximation ();
00539
00540 alpha=Tp->alpha;
00541 err=Tp->err;
00542 ini=Tp->nii;
00543
00544
00545 Tp->time=Tp->timecont.starttime ();
00546
00547 end_time = Tp->timecont.endtime ();
00548
00549 dt = Tp->timecont.initialtimeincr ();
00550
00551
00552
00553
00554 i=0;
00555
00556 Tm->initmaterialmodels();
00557 print_initt(-1, "wt");
00558 print_stept(0,i,Tp->time,NULL);
00559 do{
00560
00561 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00562
00563 Tm->updateipval ();
00564 i++;
00565
00566
00567 capacity_matrix (0);
00568
00569
00570 conductivity_matrix (0);
00571
00572
00573
00574 for (j=0;j<n;j++){
00575 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00576 }
00577
00578
00579 Kmat->gmxv (d,p);
00580
00581
00582
00583
00584
00585 Kmat->scalgm (dt*alpha);
00586 Kmat->addgm (1.0,*Cmat);
00587
00588
00589 trfel_right_hand_side (0,rhs,n);
00590
00591
00592 for (j=0;j<n;j++){
00593 fb[j] = rhs[j] - p[j];
00594 }
00595
00596
00597
00598 Tp->ssle->solve_system (Gtt,Kmat,tdlhs,fb,Outt);
00599
00600
00601 for (j=0;j<n;j++){
00602 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00603 }
00604
00605
00606 approximation ();
00607
00608
00609 internal_fluxes (fi,n);
00610
00611 stop = norm_computation (fi,rhs,err,1,1);
00612
00613 if (stop==0){
00614
00615 for (j=0;j<ini;j++){
00616
00617
00618 Tp->ssle->solve_system (Gtt,Kmat,p,fi,Outt);
00619
00620 for (k=0;k<n;k++){
00621 tdlhs[k]+=p[k];
00622 lhs[k]+=alpha*dt*p[k];
00623 }
00624
00625 approximation ();
00626
00627 internal_fluxes (fi,n);
00628
00629 if (Mesprt==1) fprintf (stdout,"\n iteration %ld",j);
00630 stop = norm_computation (fi,rhs,err,1,1);
00631
00632 if (stop==1){
00633 break;
00634 }
00635 }
00636 }
00637
00638
00639
00640 Tp->time = Tp->timecont.newtime ();
00641 dt = Tp->timecont.actualbacktimeincr ();
00642
00643 print_stept(0,i,Tp->time,NULL);
00644 print_flusht();
00645
00646 }while(Tp->time<end_time);
00647
00648 delete [] fi; delete [] fb; delete [] p;
00649 delete [] arhs; delete [] d;
00650
00651 print_closet();
00652 }
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663 void nonlinear_nonstat_solv_dform ()
00664 {
00665 long i,j,k,n,ini,stop;
00666 double dt,end_time,alpha,err;
00667 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*prhs,*plhs,*ptdlhs;
00668
00669
00670 n=Ndoft;
00671
00672
00673 lhs = Lsrst->give_lhs (0);
00674
00675 tdlhs = Lsrst->give_tdlhs (0);
00676
00677 lhsi = Lsrst->lhsi;
00678
00679 rhs = Lsrst->give_rhs (0);
00680
00681
00682 d = new double [n];
00683 prhs = new double [n];
00684 plhs = new double [n];
00685 ptdlhs = new double [n];
00686 p = new double [n];
00687 fb = new double [n];
00688 fi = new double [n];
00689
00690
00691 nullv (lhs,n);
00692 nullv (plhs,n);
00693 nullv (tdlhs,n);
00694 nullv (ptdlhs,n);
00695 nullv (rhs,n);
00696 nullv (prhs,n);
00697 nullv (d,n);
00698 nullv (p,n);
00699 nullv (fb,n);
00700 nullv (fi,n);
00701
00702
00703 approximation ();
00704
00705 alpha=Tp->alpha;
00706 err=Tp->err;
00707 ini=Tp->nii;
00708
00709
00710 Tp->time=Tp->timecont.starttime ();
00711
00712 end_time = Tp->timecont.endtime ();
00713
00714 dt = Tp->timecont.initialtimeincr ();
00715
00716
00717
00718
00719 i=0;
00720 Tm->initmaterialmodels();
00721 print_initt(-1, "wt");
00722 print_stept(0,i,Tp->time,NULL);
00723 do{
00724
00725 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00726
00727 Tm->updateipval ();
00728 i++;
00729
00730
00731 capacity_matrix (0);
00732
00733
00734 conductivity_matrix (0);
00735
00736
00737
00738
00739
00740
00741
00742 Cmat->gmxv (lhs,p);
00743
00744
00745
00746 Kmat->scalgm (dt*alpha);
00747 Kmat->addgm (1.0,*Cmat);
00748
00749
00750
00751 trfel_right_hand_side (0,rhs,n);
00752
00753 fprintf (Outt,"\n\n\n rhs %lf",rhs[0]);
00754
00755
00756 for (j=0;j<n;j++){
00757 fb[j] = rhs[j]*alpha*dt + p[j];
00758 }
00759
00760 copyv (lhs,plhs,n);
00761
00762 fprintf (Outt,"\n fb %lf\n",fb[0]);
00763
00764
00765
00766 Tp->ssle->solve_system (Gtt,Kmat,lhs,fb,Outt);
00767
00768 if (alpha<Tp->zero){
00769 for (j=0;j<n;j++){
00770 tdlhs[j]=ptdlhs[j];
00771 }
00772 }
00773 else{
00774 for (j=0;j<n;j++){
00775 tdlhs[j]=(lhs[j]-plhs[j])/dt/alpha;
00776 }
00777 }
00778
00779
00780 approximation ();
00781
00782
00783 internal_fluxes (fi,n);
00784
00785 stop = norm_computation (fi,rhs,err,1,1);
00786
00787 for (j=0;j<n;j++){
00788
00789
00790 fb[j]=0.0-fi[j];
00791 }
00792
00793 if (stop==0){
00794
00795 for (j=0;j<ini;j++){
00796
00797
00798
00799
00800 Tp->ssle->solve_system (Gtt,Kmat,p,fb,Outt);
00801
00802 for (k=0;k<n;k++){
00803 lhs[k]-=p[k];
00804 }
00805
00806 if (alpha<Tp->zero){
00807 for (k=0;k<n;k++){
00808 tdlhs[j]=ptdlhs[j];
00809 }
00810 }
00811 else{
00812 for (k=0;k<n;k++){
00813 tdlhs[k]=(lhs[k]-plhs[k])/dt/alpha;
00814 }
00815 }
00816
00817 approximation ();
00818
00819 internal_fluxes (fi,n);
00820
00821 for (k=0;k<n;k++){
00822
00823 fb[k]=rhs[k]-fi[k];
00824 }
00825
00826 if (Mesprt==1) fprintf (stdout,"\n iteration %ld",j);
00827 stop = norm_computation (fi,rhs,err,1,1);
00828
00829
00830
00831 if (stop==1){
00832 break;
00833 }
00834 }
00835 }
00836
00837
00838 copyv (rhs,prhs,n);
00839 copyv (tdlhs,ptdlhs,n);
00840
00841
00842 Tp->time = Tp->timecont.newtime ();
00843 dt = Tp->timecont.actualbacktimeincr ();
00844
00845 print_stept(0,i,Tp->time,NULL);
00846 print_flusht();
00847
00848 }while(Tp->time<end_time);
00849
00850 delete [] fi; delete [] fb; delete [] p;
00851 delete [] prhs; delete [] tdlhs; delete [] d;
00852
00853 print_closet();
00854 }
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864 void nonlinear_nonstat_solv_pokus ()
00865 {
00866 long i,j,k,n,ini;
00867 double dt,end_time,alpha,norfb,err;
00868 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*arhs;
00869 double norfb_last,*lhs_last;
00870
00871 n=Ndoft;
00872
00873 lhs = Lsrst->give_lhs (0);
00874 tdlhs = Lsrst->give_tdlhs (0);
00875 lhsi = Lsrst->lhsi;
00876 rhs = Lsrst->give_rhs (0);
00877
00878 d = new double [n];
00879 lhs_last = new double [n];
00880 arhs = new double [n];
00881 p = new double [n];
00882 fb = new double [n];
00883 fi = new double [n];
00884
00885
00886 nullv (lhs,n);
00887 nullv (lhs_last,n);
00888 nullv (tdlhs,n);
00889 nullv (rhs,n);
00890 nullv (arhs,n);
00891 nullv (d,n);
00892 nullv (p,n);
00893 nullv (fb,n);
00894 nullv (fi,n);
00895
00896
00897 approximation ();
00898
00899 alpha=Tp->alpha;
00900 err=Tp->err;
00901 ini=Tp->nii;
00902
00903
00904 Tp->time=Tp->timecont.starttime ();
00905
00906 end_time = Tp->timecont.endtime ();
00907
00908 dt = Tp->timecont.initialtimeincr ();
00909
00910
00911
00912
00913 i=0;
00914
00915 Tm->initmaterialmodels();
00916 print_initt(-1, "wt");
00917 print_stept(0,i,Tp->time,NULL);
00918 do{
00919
00920 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
00921
00922 Tm->updateipval ();
00923
00924
00925 capacity_matrix (0);
00926
00927
00928 conductivity_matrix (0);
00929
00930 for (j=0;j<n;j++){
00931 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00932 }
00933
00934
00935 trfel_right_hand_side (0,rhs,n);
00936
00937
00938 Kmat->gmxv (d,p);
00939
00940
00941
00942 Kmat->scalgm (dt*alpha);
00943 Kmat->addgm (1.0,*Cmat);
00944
00945
00946 for (j=0;j<n;j++){
00947 fb[j] = rhs[j] - p[j];
00948 }
00949
00950
00951 Tp->ssle->solve_system (Gtt,Kmat,tdlhs,fb,Outt);
00952
00953
00954 for (j=0;j<n;j++){
00955 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00956 }
00957
00958 approximation ();
00959 internal_fluxes (fi,n);
00960
00961
00962 for (j=0;j<n;j++){
00963 fb[j]=fi[j];
00964
00965 }
00966
00967 norfb = ss (fb,fb,n);
00968
00969
00970 for (k=0;k<n;k++){
00971 lhs_last[k] = lhs[k];
00972 }
00973 norfb_last = norfb;
00974
00975
00976 if (Mesprt==1) fprintf (stdout,"\n%e %e",norfb,err);
00977
00978
00979 if (norfb<err){
00980
00981 Tp->time = Tp->timecont.newtime ();
00982 dt = Tp->timecont.actualbacktimeincr ();
00983
00984 i++;
00985 print_stept(0,i,Tp->time,NULL);
00986 print_flusht();
00987 continue;
00988
00989 }
00990
00991 for (j=0;j<ini;j++){
00992
00993
00994
00995
00996
00997
00998 capacity_matrix (0);
00999
01000
01001 conductivity_matrix (0);
01002
01003 Kmat->gmxv (d,p);
01004
01005
01006
01007 Kmat->scalgm (dt*alpha);
01008 Kmat->addgm (1.0,*Cmat);
01009
01010
01011
01012 Tp->ssle->solve_system (Gtt,Kmat,p,fb,Outt);
01013
01014 for (k=0;k<n;k++){
01015 tdlhs[k]+=p[k];
01016 lhs[k]+=alpha*dt*p[k];
01017 }
01018
01019 approximation ();
01020 internal_fluxes (fi,n);
01021
01022
01023 for (k=0;k<n;k++){
01024 fb[k]= fi[k];
01025 }
01026
01027 norfb = ss (fb,fb,n);
01028
01029 if (Mesprt==1) fprintf (stdout,"\n iteration %ld error %e",j,norfb);
01030
01031 approximation ();
01032
01033 if (norfb<err){
01034 break;
01035 }
01036
01037
01038
01039
01040
01041 if (norfb>norfb_last){
01042 for (k=0;k<n;k++){
01043 lhs[k] = lhs_last[k];
01044 }
01045 approximation ();
01046 break;
01047 }
01048 for (k=0;k<n;k++){
01049 lhs_last[k] = lhs[k];
01050 }
01051
01052 norfb_last = norfb;
01053
01054
01055 }
01056
01057
01058
01059 Tp->time = Tp->timecont.newtime ();
01060 dt = Tp->timecont.actualbacktimeincr ();
01061 i++;
01062
01063 print_stept(0,i,Tp->time,NULL);
01064 print_flusht();
01065
01066 }while(Tp->time<end_time);
01067
01068 delete [] fi; delete [] fb; delete [] p;
01069 delete [] arhs; delete [] d; delete [] lhs_last;
01070
01071 print_closet();
01072 }
01073
01074
01075
01076
01077
01078
01079
01080 void nonlinear_nonstat_solv_new ()
01081 {
01082 long i,j,k,n,ini,stop;
01083 double dt,end_time,alpha,norfb,norrhs,err;
01084 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*arhs;
01085 double *lhs_last;
01086
01087 n=Ndoft;
01088
01089 lhs = Lsrst->give_lhs (0);
01090 tdlhs = Lsrst->give_tdlhs (0);
01091 lhsi = Lsrst->lhsi;
01092 rhs = Lsrst->give_rhs (0);
01093
01094 d = new double [n];
01095 lhs_last = new double [n];
01096 arhs = new double [n];
01097 p = new double [n];
01098 fb = new double [n];
01099 fi = new double [n];
01100
01101
01102 nullv (lhs,n);
01103 nullv (lhs_last,n);
01104 nullv (tdlhs,n);
01105 nullv (rhs,n);
01106 nullv (arhs,n);
01107 nullv (d,n);
01108 nullv (p,n);
01109 nullv (fb,n);
01110 nullv (fi,n);
01111
01112
01113 approximation ();
01114
01115 alpha=Tp->alpha;
01116 err=Tp->err;
01117 ini=Tp->nii;
01118
01119
01120 Tp->time=Tp->timecont.starttime ();
01121
01122 end_time = Tp->timecont.endtime ();
01123
01124 dt = Tp->timecont.initialtimeincr ();
01125
01126
01127
01128
01129 i=0;
01130
01131 Tm->initmaterialmodels();
01132 print_initt(-1, "wt");
01133 print_stept(0,i,Tp->time,NULL);
01134 do{
01135
01136 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
01137
01138 Tm->updateipval ();
01139 i++;
01140
01141
01142 capacity_matrix (0);
01143
01144
01145
01146
01147
01148
01149 conductivity_matrix (0);
01150
01151 for (j=0;j<n;j++){
01152 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
01153 }
01154
01155
01156 trfel_right_hand_side (0,rhs,n);
01157
01158
01159 norrhs = ss (rhs,rhs,n);
01160
01161
01162 Kmat->gmxv (d,p);
01163
01164
01165
01166 Kmat->scalgm (dt*alpha);
01167 Kmat->addgm (1.0,*Cmat);
01168
01169
01170 for (j=0;j<n;j++){
01171 fb[j] = rhs[j] - p[j];
01172 }
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186 fprintf (Outt,"\n");
01187
01188
01189
01190
01191
01192 Tp->ssle->solve_system (Gtt,Kmat,tdlhs,fb,Outt);
01193
01194
01195
01196
01197 for (j=0;j<n;j++){
01198
01199
01200 lhs[j] = d[j] + alpha*dt*tdlhs[j];
01201
01202 }
01203
01204
01205
01206 for (j=0;j<ini;j++){
01207
01208
01209 solution_correction ();
01210 approximation ();
01211
01212 fprintf (Outt,"\n\n krok %ld",i);
01213 fprintf (Outt,"\n %le %le %le %le",Tm->ip[0].av[0],Tm->ip[1].av[0],Tm->ip[4].av[0],Tm->ip[5].av[0]);
01214
01215
01216
01217
01218
01219 if (Tp->trsolv == fullnewtont){
01220
01221 capacity_matrix (0);
01222
01223
01224 conductivity_matrix (0);
01225
01226
01227 Kmat->gmxv (d,p);
01228
01229
01230
01231 Kmat->scalgm (dt*alpha);
01232 Kmat->addgm (1.0,*Cmat);
01233
01234 }
01235
01236
01237
01238 solution_correction ();
01239 approximation ();
01240
01241
01242 if (Tp->trestype==lrhst){
01243 Kmat->gmxv (tdlhs,fi);
01244
01245 for (k=0;k<n;k++){
01246 fb[k] = rhs[k] - p[k] - fi[k];
01247 }
01248 }
01249 if (Tp->trestype==fluxest){
01250 internal_fluxes (fi,n);
01251
01252 for (k=0;k<n;k++){
01253 fb[k]=fi[k];
01254
01255 }
01256 }
01257
01258
01259 norfb = ss (fb,fb,n);
01260
01261 stop = norm_computation (fb,rhs,err,i,j);
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281 if (stop==Tp->ntm)
01282 break;
01283
01284
01285 Tp->ssle->solve_system (Gtt,Kmat,p,fb,Outt);
01286
01287
01288 for (k=0;k<n;k++){
01289
01290
01291 tdlhs[k]+=p[k];
01292 lhs[k]+=alpha*dt*p[k];
01293
01294 }
01295
01296
01297 }
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318 Tp->time = Tp->timecont.newtime ();
01319 dt = Tp->timecont.actualbacktimeincr ();
01320
01321 print_stept(0,i,Tp->time,NULL);
01322 print_flusht();
01323
01324 }while(Tp->time<end_time);
01325
01326 delete [] fi; delete [] fb; delete [] p;
01327 delete [] arhs; delete [] d; delete [] lhs_last;
01328
01329 print_closet();
01330 }
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342 void nonlinear_nonstat_solv_linesearch ()
01343 {
01344 long i,j,k,l,n,ini,nils;
01345 double dt,end_time,alpha,norfb,err,eta,deta,a,errls,nom;
01346 double *d,*fb,*fi,*p,*pp,*lhs,*lhsi,*tdlhs,*rhs;
01347
01348 nils=10;
01349 errls=1.0e-6;
01350
01351 n=Ndoft;
01352
01353 lhs = Lsrst->give_lhs (0);
01354 tdlhs = Lsrst->give_tdlhs (0);
01355 lhsi = Lsrst->lhsi;
01356 rhs = Lsrst->give_rhs (0);
01357
01358 d = new double [n];
01359 p = new double [n];
01360 pp = new double [n];
01361 fb = new double [n];
01362 fi = new double [n];
01363
01364
01365 nullv (lhs,n);
01366 nullv (tdlhs,n);
01367 nullv (rhs,n);
01368 nullv (d,n);
01369 nullv (p,n);
01370 nullv (fb,n);
01371 nullv (fi,n);
01372
01373
01374
01375 approximation ();
01376
01377 alpha=Tp->alpha;
01378 err=Tp->err;
01379 ini=Tp->nii;
01380
01381
01382 Tp->time=Tp->timecont.starttime ();
01383
01384 end_time = Tp->timecont.endtime ();
01385
01386 dt = Tp->timecont.initialtimeincr ();
01387
01388
01389
01390
01391 i=0;
01392 Tm->initmaterialmodels();
01393 print_initt(-1, "wt");
01394 print_stept(0,i,Tp->time,NULL);
01395 do{
01396
01397 if (Mesprt != 0) fprintf (stdout,"\n\n iteration number %ld, time %f, time step = %f\n",i,Tp->time,dt);
01398
01399 Tm->updateipval ();
01400
01401
01402 capacity_matrix (0);
01403
01404
01405 conductivity_matrix (0);
01406
01407
01408 for (j=0;j<n;j++){
01409 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
01410 }
01411
01412
01413 trfel_right_hand_side (0,rhs,n);
01414
01415
01416 Kmat->gmxv (d,p);
01417
01418
01419
01420 Kmat->scalgm (dt*alpha);
01421 Kmat->addgm (1.0,*Cmat);
01422
01423
01424 for (j=0;j<n;j++){
01425 fb[j] = rhs[j] - p[j];
01426 }
01427
01428
01429
01430 Tp->ssle->solve_system (Gtt,Kmat,p,fb,Outt);
01431
01432
01433 for (j=0;j<n;j++){
01434 tdlhs[j]=p[j];
01435 lhs[j] = d[j] + alpha*dt*tdlhs[j];
01436 }
01437
01438 approximation ();
01439
01440
01441 internal_fluxes (fb,n);
01442
01443
01444 norfb = ss (fb,fb,n);
01445
01446 if (Mesprt==1) fprintf (stdout,"\n%e %e",norfb,err);
01447
01448
01449 if (norfb<err){
01450
01451 Tp->time = Tp->timecont.newtime ();
01452 dt = Tp->timecont.actualbacktimeincr ();
01453
01454 i++;
01455 print_stept(0,i,Tp->time,NULL);
01456 print_flusht();
01457 continue;
01458
01459 }
01460
01461
01462 for (j=0;j<ini;j++){
01463
01464
01465
01466 capacity_matrix (0);
01467
01468 conductivity_matrix (0);
01469 Kmat->scalgm (dt*alpha);
01470 Kmat->addgm (1.0,*Cmat);
01471
01472
01473
01474 Kmat->gmxv (p,pp);
01475 a = ss (pp,p,n);
01476 eta =1.0;
01477
01478 fprintf (stdout,"\n jsme ve vnitrni smycce");
01479
01480 for (k=0;k<nils;k++){
01481 nom = ss (fb,p,n);
01482 deta = nom/a;
01483 eta-=deta;
01484
01485
01486 for (l=0;l<n;l++){
01487 tdlhs[l]=eta*p[l];
01488 lhs[l]=d[l]+alpha*dt*tdlhs[l];
01489 }
01490 approximation ();
01491 internal_fluxes (fb,n);
01492
01493 norfb = ss (fb,fb,n);
01494 fprintf (stdout,"\n jsme v line search, nom=%le deta %le eta %le norfb %le",nom,deta,eta,norfb);
01495
01496 if (fabs(nom)<errls)
01497 break;
01498
01499 }
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513 norfb = ss (fb,fb,n);
01514
01515 if (Mesprt==1) fprintf (stdout,"\n iteration %ld error %e",j,norfb);
01516
01517 if (norfb<err){
01518 break;
01519 }
01520
01521 Tp->ssle->solve_system (Gtt,Kmat,p,fb,Outt);
01522
01523 }
01524
01525
01526 Tp->time = Tp->timecont.newtime ();
01527 dt = Tp->timecont.actualbacktimeincr ();
01528 i++;
01529
01530 print_stept(0,i,Tp->time,NULL);
01531 print_flusht();
01532
01533 }while(Tp->time<end_time);
01534
01535 delete [] fi; delete [] fb; delete [] p;
01536 delete [] d;
01537
01538 print_closet();
01539 }
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568 long norm_computation (double *f,double *rhs,double err,long nt,long sc)
01569 {
01570 long i,j,dof,stop,pstop;
01571 double norf,norrhs;
01572
01573 if (nt==1){
01574
01575
01576
01577 norf=0.0;
01578 for (i=0;i<Ndoft;i++){
01579 norf+=f[i]*f[i];
01580 }
01581 norf=sqrt(norf);
01582
01583 if (sc==0){
01584
01585 if (norf<err)
01586 stop=1;
01587 else
01588 stop=0;
01589
01590 if (Mesprt==1) fprintf (stdout,"\n norm of residuum %e, required error %e",norf,err);
01591 }
01592
01593 if (sc==1){
01594
01595 norrhs=0.0;
01596 for (i=0;i<Ndoft;i++){
01597 norrhs+=rhs[i]*rhs[i];
01598 }
01599 norrhs=sqrt(norrhs);
01600
01601 if (norrhs>0.0){
01602 if (norf/norrhs<err)
01603 stop=1;
01604 else
01605 stop=0;
01606
01607 if (Mesprt==1) fprintf (stdout,"\n norm of residuum %e, norm of rhs %e, error %e, required error %e",norf,norrhs,norf/norrhs,err);
01608 }
01609 else{
01610 if (norf<err)
01611 stop=1;
01612 else
01613 stop=0;
01614
01615 if (Mesprt==1) fprintf (stdout,"\n norm of residuum %e, norm of rhs %e, required error %e",norf,norrhs,err);
01616 }
01617 }
01618 }
01619
01620
01621 if (nt==2){
01622
01623
01624
01625
01626
01627 pstop=0;
01628 for (i=0;i<Tp->ntm;i++){
01629 norf=0.0;
01630 norrhs=0.0;
01631 for (j=0;j<Tt->nn;j++){
01632 dof = Tt->give_dof (j,i);
01633 if (dof>0){
01634 dof--;
01635
01636 norf+=(f[dof]-rhs[dof])*(f[dof]-rhs[dof]);
01637 norrhs+=rhs[dof]*rhs[dof];
01638 }
01639 norf=sqrt(norf);
01640 norrhs=sqrt(norrhs);
01641 }
01642
01643 if (sc==0){
01644
01645 if (norf<err)
01646 pstop++;
01647
01648 if (Mesprt==1) fprintf (stdout,"\n phase %ld, norm of residuum %e, required error %e",i+1,norf,err);
01649 }
01650
01651 if (sc==1){
01652
01653 if (norrhs>0.0){
01654 if (norf/norrhs<err)
01655 pstop++;
01656
01657 if (Mesprt==1) fprintf (stdout,"\n phase %ld, norm of residuum %e, norm of rhs %e, error %e, required error %e",i+1,norf,norrhs,norf/norrhs,err);
01658 }
01659 else{
01660 if (norf<err)
01661 pstop++;
01662
01663 if (Mesprt==1) fprintf (stdout,"\n phase %ld, norm of residuum %e, norm of rhs %e, required error %e",i+1,norf,norrhs,err);
01664 }
01665 }
01666 }
01667
01668 if (pstop==Tp->ntm)
01669 stop=1;
01670 else
01671 stop=0;
01672 }
01673
01674 return stop;
01675 }
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705 long norm_computation_vec (double *f,double *rhs,double *err,double *thresh,long nt,long sc)
01706 {
01707 long i,j,dof,stop,pstop;
01708 double norf,norrhs;
01709
01710 if (nt==1){
01711
01712
01713
01714 norf=0.0;
01715 for (i=0;i<Ndoft;i++){
01716 norf+=f[i]*f[i];
01717 }
01718 norf=sqrt(norf);
01719
01720 if (sc==0){
01721
01722 if (norf<err[0])
01723 stop=1;
01724 else
01725 stop=0;
01726
01727 if (Mesprt==1) fprintf (stdout,"\n norm of residuum %e, required error %e",norf,err[0]);
01728 }
01729
01730 if (sc==1){
01731
01732 norrhs=0.0;
01733 for (i=0;i<Ndoft;i++){
01734 norrhs+=rhs[i]*rhs[i];
01735 }
01736 norrhs=sqrt(norrhs);
01737
01738 if (norrhs>thresh[0]){
01739 if (norf/norrhs<err[0])
01740 stop=1;
01741 else
01742 stop=0;
01743
01744 if (Mesprt==1) fprintf (stdout,"\n norm of residuum %e, norm of rhs %e, error %e, required error %e",norf,norrhs,norf/norrhs,err[0]);
01745 }
01746 else{
01747 if (norf<err[0])
01748 stop=1;
01749 else
01750 stop=0;
01751
01752 if (Mesprt==1) fprintf (stdout,"\n norm of residuum %e, norm of rhs %e, required error %e",norf,norrhs,err[0]);
01753 }
01754 }
01755 }
01756
01757
01758 if (nt==2){
01759
01760
01761
01762
01763
01764 pstop=0;
01765 for (i=0;i<Tp->ntm;i++){
01766 norf=0.0;
01767 norrhs=0.0;
01768 for (j=0;j<Tt->nn;j++){
01769 dof = Tt->give_dof (j,i);
01770 if (dof>0){
01771 dof--;
01772 norf+=f[dof]*f[dof];
01773
01774 norrhs+=rhs[dof]*rhs[dof];
01775 }
01776 }
01777 norf=sqrt(norf);
01778 norrhs=sqrt(norrhs);
01779
01780 if (sc==0){
01781
01782 if (norf<err[i])
01783 pstop++;
01784
01785 if (Mesprt==1) fprintf (stdout,"\n phase %ld, norm of residuum %e, required error %e",i+1,norf,err[i]);
01786 }
01787
01788 if (sc==1){
01789
01790 if (norrhs>thresh[i]){
01791 if (norf/norrhs<err[i])
01792 pstop++;
01793
01794 if (Mesprt==1) fprintf (stdout,"\n phase %ld, norm of residuum %e, norm of rhs %e, error %e, required error %e",i+1,norf,norrhs,norf/norrhs,err[i]);
01795 }
01796 else{
01797 if (norf<err[i])
01798 pstop++;
01799
01800 if (Mesprt==1) fprintf (stdout,"\n phase %ld, norm of residuum %e, norm of rhs %e, required error %e",i+1,norf,norrhs,err[i]);
01801 }
01802 }
01803
01804 }
01805
01806 if (pstop==Tp->ntm)
01807 stop=1;
01808 else
01809 stop=0;
01810 }
01811
01812 return stop;
01813 }
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868 void nonlinear_nonstat_solv_dform_dneska ()
01869 {
01870 long i,j,k,n,ini;
01871 double dt,end_time,alpha,err,norm;
01872 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*plhs,*pplhs,*rhs;
01873
01874
01875 n=Ndoft;
01876
01877
01878 lhs = Lsrst->give_lhs (0);
01879
01880 tdlhs = Lsrst->give_tdlhs (0);
01881
01882 lhsi = Lsrst->lhsi;
01883
01884 plhs = new double [n];
01885
01886 pplhs = new double [n];
01887
01888 rhs = Lsrst->give_rhs (0);
01889
01890
01891
01892 d = new double [n];
01893
01894 p = new double [n];
01895
01896 fb = new double [n];
01897
01898 fi = new double [n];
01899
01900
01901 nullv (lhs,n);
01902 nullv (plhs,n);
01903 nullv (rhs,n);
01904 nullv (d,n);
01905 nullv (fb,n);
01906
01907 for (i=0;i<n;i++){
01908 plhs[i]=lhsi[i];
01909 pplhs[i]=lhsi[i];
01910 }
01911
01912
01913
01914
01915
01916
01917
01918 alpha=Tp->alpha;
01919 err=Tp->err;
01920 ini=Tp->nii;
01921
01922
01923 Tp->time=Tp->timecont.starttime ();
01924
01925 end_time = Tp->timecont.endtime ();
01926
01927 dt = Tp->timecont.initialtimeincr ();
01928
01929
01930
01931
01932 i=0;
01933
01934
01935 approximation ();
01936 if (Tp->nvs==1 && Tp->pnvs==1)
01937 actual_previous_change ();
01938
01939 Tm->initmaterialmodels();
01940 print_initt(-1, "wt");
01941 print_stept(0,i,Tp->time,NULL);
01942
01943 do{
01944
01945 if (Mesprt != 0) fprintf (stdout,"\n\n increment number %ld, time %f, time step = %f\n",i,Tp->time,dt);
01946
01947 Tm->updateipval ();
01948 i++;
01949
01950
01951 capacity_matrix (0);
01952
01953
01954 conductivity_matrix (0);
01955
01956
01957
01958
01959
01960
01961
01962
01963 nullv (fb,n);
01964 if (Tp->tprob == discont_nonstat_problem || Tp->tprob == discont_nonlin_nonstat_problem){
01965
01966 }
01967
01968
01969 for (j=0;j<n;j++){
01970 d[j]=lhs[j]+(1.0-alpha)*dt*tdlhs[j];
01971 }
01972
01973
01974 Cmat->gmxv (d,p);
01975 cmulv(1.0/alpha/dt,p,n);
01976
01977
01978
01979
01980 Cmat->scalgm (1.0/dt/alpha);
01981 Kmat->addgm (1.0,*Cmat);
01982
01983
01984
01985 trfel_right_hand_side (0,rhs,n);
01986
01987
01988 for (j=0;j<n;j++){
01989 fb[j] += rhs[j] + p[j];
01990 }
01991
01992
01993 Tp->ssle->solve_system (Gtt,Kmat,lhs,fb,Outt);
01994
01995 for (j=0;j<n;j++){
01996 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
01997 }
01998
01999
02000 solution_correction ();
02001
02002 approximation ();
02003
02004
02005
02006
02007 capacity_matrix (0);
02008
02009
02010 conductivity_matrix (0);
02011
02012 Cmat->gmxv (tdlhs,p);
02013 Kmat->gmxv (lhs,fi);
02014
02015 nullv (fb,n);
02016 if (Tp->tprob == discont_nonstat_problem || Tp->tprob == discont_nonlin_nonstat_problem){
02017
02018 }
02019
02020 norm=0.0;
02021 for (j=0;j<n;j++){
02022 fi[j] = rhs[j] + fb[j] - fi[j] - p[j];
02023 norm+=fi[j]*fi[j];
02024 }
02025
02026
02027
02028
02029
02030 if (norm>err){
02031
02032 for (j=0;j<ini;j++){
02033
02034
02035
02036
02037 Cmat->scalgm (1.0/dt/alpha);
02038 Kmat->addgm (1.0,*Cmat);
02039
02040 Tp->ssle->solve_system (Gtt,Kmat,p,fi,Outt);
02041
02042 for (k=0;k<n;k++){
02043 lhs[k]+=p[k];
02044 }
02045
02046 for (k=0;k<n;k++){
02047 tdlhs[k]=(lhs[k]-d[k])/dt/alpha;
02048 }
02049
02050
02051 solution_correction ();
02052 approximation ();
02053
02054
02055
02056 capacity_matrix (0);
02057
02058
02059 conductivity_matrix (0);
02060
02061 Cmat->gmxv (tdlhs,p);
02062 Kmat->gmxv (lhs,fi);
02063
02064 nullv (fb,n);
02065 if (Tp->tprob == discont_nonstat_problem || Tp->tprob == discont_nonlin_nonstat_problem){
02066
02067 }
02068
02069 norm=0.0;
02070 for (k=0;k<n;k++){
02071 fi[k] = rhs[k] + fb[k] - fi[k] - p[k];
02072 norm+=fi[k]*fi[k];
02073 }
02074
02075
02076
02077
02078
02079
02080 if (norm<err){
02081 break;
02082 }
02083 }
02084 }
02085
02086
02087 compute_cycles (plhs,pplhs);
02088
02089 if (i>1){
02090 for (j=0;j<n;j++){
02091 pplhs[j]=plhs[j];
02092 }
02093 }
02094 for (j=0;j<n;j++){
02095 plhs[j]=lhs[j];
02096 }
02097
02098
02099
02100
02101
02102
02103
02104 Tp->time = Tp->timecont.newtime ();
02105 dt = Tp->timecont.actualbacktimeincr ();
02106
02107 print_stept(0,i,Tp->time,NULL);
02108 print_flusht();
02109
02110 }while(Tp->time<end_time);
02111
02112 fprintf(Outt, "\n\n\n pocty cyklu \n");
02113 for (i = 0; i <Tt->nn;i++){
02114 fprintf(Outt, "\n node %6ld %ld",i,Tt->ncycl[i]);
02115 }
02116
02117 delete [] fb; delete [] p;
02118 delete [] d;
02119
02120 print_closet();
02121 }
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170
02171
02172
02173
02174 void nonlinear_nonstat_solv_fnr_dform_old ()
02175 {
02176 long i,j,k,n,ini,ani,nsts,lcid;
02177 double dt,dtdef,dtmin,end_time,time,alpha,err,norres,norrhs;
02178 double *d,*f,*p,*v,*lhs,*lhsi,*tdlhs,*rhs,*lhsb,*tdlhsb;
02179
02180
02181 n=Ndoft;
02182
02183 lcid=0;
02184
02185
02186 lhs = Lsrst->give_lhs (0);
02187
02188 tdlhs = Lsrst->give_tdlhs (0);
02189
02190 lhsi = Lsrst->lhsi;
02191
02192 rhs = Lsrst->give_rhs (0);
02193
02194
02195
02196 d = new double [n];
02197
02198 p = new double [n];
02199
02200 f = new double [n];
02201
02202 v = new double [n];
02203
02204 lhsb = new double [n];
02205
02206 tdlhsb = new double [n];
02207
02208
02209
02210 nullv (lhs,n);
02211 nullv (rhs,n);
02212 nullv (tdlhs,n);
02213 nullv (d,n);
02214 nullv (p,n);
02215 nullv (f,n);
02216 nullv (v,n);
02217
02218
02219 approximation ();
02220
02221
02222 alpha=Tp->alpha;
02223
02224 err=Tp->err;
02225
02226 ini=Tp->nii;
02227
02228 dtmin=Tp->timecont.dtmin;
02229
02230
02231 Tp->time=Tp->timecont.starttime ();
02232 time=Tp->time;
02233
02234 end_time = Tp->timecont.endtime ();
02235
02236 dt = Tp->timecont.initialtimeincr ();
02237
02238
02239
02240
02241
02242 i=0;
02243
02244 nsts=0;
02245
02246 Tm->initmaterialmodels();
02247 print_initt(-1, "wt");
02248 print_stept(0,i,Tp->time,NULL);
02249 do{
02250
02251 if (Mesprt != 0) fprintf (stdout,"\n\n time %f, time step = %f, increment number %ld",Tp->time,dt,i);
02252
02253 Tm->updateipval ();
02254 i++;
02255
02256
02257 for (j=0;j<n;j++){
02258 lhsb[j]=lhs[j];
02259 tdlhsb[j]=tdlhs[j];
02260 }
02261
02262
02263 capacity_matrix (lcid);
02264
02265
02266 conductivity_matrix (lcid);
02267
02268
02269 for (j=0;j<n;j++){
02270 d[j] = lhs[j]+dt*(1.0-alpha)*tdlhs[j];
02271 }
02272
02273
02274 Cmat->gmxv (d,p);
02275
02276
02277
02278 Kmat->scalgm (dt*alpha);
02279 Kmat->addgm (1.0,*Cmat);
02280
02281
02282
02283 trfel_right_hand_side (lcid,f,n);
02284
02285
02286 norrhs=0.0;
02287 for (j=0;j<n;j++){
02288 rhs[j] = f[j]*alpha*dt + p[j];
02289 norrhs+=f[j]*f[j];
02290 }
02291 norrhs=sqrt(norrhs);
02292
02293
02294
02295
02296 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
02297
02298
02299 for (j=0;j<n;j++){
02300 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
02301 }
02302
02303
02304 approximation ();
02305
02306
02307
02308 for (j=0;j<ini;j++){
02309
02310
02311 capacity_matrix (lcid);
02312
02313
02314 conductivity_matrix (lcid);
02315
02316
02317
02318
02319 Kmat->scalgm (dt*alpha);
02320 Kmat->addgm (1.0,*Cmat);
02321
02322
02323 Kmat->gmxv (lhs,v);
02324
02325
02326 Cmat->gmxv (d,p);
02327
02328
02329 norres=0.0;
02330 for (k=0;k<n;k++){
02331 rhs[k]=f[k]*alpha*dt+p[k]-v[k];
02332 norres+=rhs[k]*rhs[k];
02333 }
02334 norres=sqrt(norres);
02335
02336 if (Mesprt != 0)
02337 fprintf (stdout,"\n nres %15.10le, nrhs %15.10le",norres,norrhs);
02338
02339 if (norrhs<Tp->zero){
02340 if (norres<err){
02341 break;
02342 }
02343 }
02344 else{
02345 if (Mesprt != 0)
02346 fprintf (stdout," nres/nrhs %15.10le",norres/norrhs);
02347 if (norres/norrhs<err){
02348 break;
02349 }
02350 }
02351
02352 if (Mesprt != 0) fprintf (stdout,"\n iteration number %ld",j);
02353
02354
02355 Tp->ssle->solve_system (Gtt,Kmat,v,rhs,Outt);
02356
02357
02358 for (k=0;k<n;k++){
02359 lhs[k]+=v[k];
02360 }
02361
02362
02363
02364 for (k=0;k<n;k++){
02365 tdlhs[k]=(lhs[k]-d[k])/dt/alpha;
02366 }
02367
02368
02369 approximation ();
02370
02371 }
02372
02373
02374 ani=j;
02375
02376 if (ani==0)
02377 nsts++;
02378
02379 if (ani==ini){
02380
02381 for (j=0;j<n;j++){
02382 lhs[j]=lhsb[j];
02383 tdlhs[j]=tdlhsb[j];
02384 }
02385
02386
02387 dt/=2.0;
02388 i--;
02389 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium");
02390 if (dt<dtmin){
02391 fprintf (stderr,"\n\n time increment is less than minimum time increment");
02392 fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
02393 }
02394 }else{
02395
02396 Tp->time = Tp->timecont.newtime ();
02397 dtdef = Tp->timecont.actualforwtimeincr ();
02398
02399 if (nsts==2){
02400 dt*=2.0;
02401 nsts=0;
02402 fprintf (stdout,"\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
02403 }
02404 if (dt>dtdef)
02405 dt=dtdef;
02406
02407 time+=dt;
02408 Tp->time=time;
02409 Tp->timecont.time=time;
02410
02411
02412 print_stept(0,i,Tp->time,NULL);
02413 print_flusht();
02414 }
02415
02416 }while(Tp->time<end_time);
02417
02418 delete [] v;
02419 delete [] f;
02420 delete [] p;
02421 delete [] d;
02422 delete [] lhsb;
02423 delete [] tdlhsb;
02424
02425 print_closet();
02426 }
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442 void nonlinear_nonstat_solv_fnr_dform ()
02443 {
02444 long i,j,k,n,ini,ani,nsts,lcid;
02445 long stop;
02446 double dt,dtmin,end_time,time,alpha,*err,norrhs,zero,*thresh;
02447 double *d,*f,*p,*v,*lhs,*lhsi,*tdlhs,*rhs,*lhsb,*tdlhsb;
02448 double *z,*diag;
02449
02450
02451 n=Ndoft;
02452
02453 lcid=0;
02454
02455 zero=1.0e-20;
02456
02457 thresh=Tp->threshrhs;
02458
02459
02460 lhs = Lsrst->give_lhs (0);
02461
02462 tdlhs = Lsrst->give_tdlhs (0);
02463
02464 lhsi = Lsrst->lhsi;
02465
02466 rhs = Lsrst->give_rhs (0);
02467
02468
02469
02470 d = new double [n];
02471
02472 p = new double [n];
02473
02474 f = new double [n];
02475
02476 v = new double [n];
02477
02478 lhsb = new double [n];
02479
02480 tdlhsb = new double [n];
02481
02482 z = new double [n];
02483 diag = new double [n];
02484
02485
02486
02487 nullv (lhs,n);
02488 nullv (rhs,n);
02489 nullv (tdlhs,n);
02490 nullv (d,n);
02491 nullv (p,n);
02492 nullv (f,n);
02493 nullv (v,n);
02494 nullv (z,n);
02495
02496
02497 approximation ();
02498
02499
02500 alpha=Tp->alpha;
02501
02502 err=Tp->errarr;
02503
02504 ini=Tp->nii;
02505
02506 dtmin=Tp->timecont.dtmin;
02507
02508
02509 Tp->time=Tp->timecont.starttime ();
02510 time=Tp->time;
02511
02512 end_time = Tp->timecont.endtime ();
02513
02514 dt = Tp->timecont.initialtimeincr ();
02515
02516
02517 Tm->initmaterialmodels();
02518
02519 approximation ();
02520
02521
02522
02523 if (Tp->nvs==1 && Tp->pnvs==1)
02524 actual_previous_nodval ();
02525
02526
02527
02528
02529
02530
02531
02532
02533 i=0;
02534
02535 nsts=0;
02536
02537 stop=0;
02538
02539
02540 Tm->initmaterialmodels();
02541 print_initt(-1, "wt");
02542 print_stept(0,i,Tp->time,NULL);
02543
02544
02545
02546 do{
02547
02548 Tp->time=Tp->timecont.newtime (dt);
02549
02550 if (Mesprt != 0) fprintf (stdout,"\n\n time %le, time step = %le, increment number %ld",Tp->time,dt,i);
02551
02552 Tm->updateipval ();
02553 i++;
02554
02555
02556 for (j=0;j<n;j++){
02557 lhsb[j]=lhs[j];
02558 tdlhsb[j]=tdlhs[j];
02559 }
02560
02561
02562 capacity_matrix (lcid);
02563
02564
02565
02566
02567
02568 conductivity_matrix (lcid);
02569
02570
02571
02572
02573
02574
02575
02576
02577
02578
02579 for (j=0;j<n;j++){
02580 d[j] = lhs[j]+dt*(1.0-alpha)*tdlhs[j];
02581 }
02582
02583
02584 Cmat->gmxv (d,p);
02585
02586
02587
02588 Kmat->scalgm (dt*alpha);
02589 Kmat->addgm (1.0,*Cmat);
02590
02591
02592
02593 trfel_right_hand_side (lcid,f,n);
02594
02595
02596 norrhs=0.0;
02597 for (j=0;j<n;j++){
02598 rhs[j] = f[j]*alpha*dt + p[j];
02599
02600
02601
02602
02603 }
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613 Kmat->diag_check (zero,rhs);
02614
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
02625
02626
02627
02628
02629
02630
02631
02632
02633 for (j=0;j<n;j++){
02634 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
02635 }
02636
02637
02638
02639
02640
02641
02642
02643 solution_correction ();
02644
02645 approximation ();
02646
02647
02648
02649 if (Tp->nvs==1 && Tp->pnvs==1)
02650 actual_previous_nodval ();
02651
02652
02653
02654 for (j=0;j<ini;j++){
02655
02656
02657 capacity_matrix (lcid);
02658
02659
02660 conductivity_matrix (lcid);
02661
02662
02663
02664
02665 Kmat->scalgm (dt*alpha);
02666 Kmat->addgm (1.0,*Cmat);
02667
02668
02669 Kmat->gmxv (lhs,v);
02670
02671
02672 Cmat->gmxv (d,p);
02673
02674
02675
02676
02677 for (k=0;k<n;k++){
02678 rhs[k]=f[k]*alpha*dt+p[k]-v[k];
02679
02680
02681 z[k]=f[k]*alpha*dt+p[k];
02682
02683
02684 }
02685
02686
02687
02688
02689
02690 stop = norm_computation_vec (rhs,z,err,thresh,2,1);
02691
02692
02693
02694 if (stop==1){
02695 break;
02696 }
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725 if (Mesprt != 0) fprintf (stdout,"\n iteration number %ld",j);
02726
02727
02728 for (k=0;k<n;k++){
02729 rhs[k]=f[k]*alpha*dt+p[k];
02730
02731 }
02732
02733
02734
02735
02736
02737
02738 Kmat->diag_check (zero,rhs);
02739
02740
02741
02742
02743
02744
02745
02746
02747 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
02748
02749
02750
02751
02752
02753
02754
02755
02756
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766
02767
02768
02769 for (k=0;k<n;k++){
02770 tdlhs[k]=(lhs[k]-d[k])/dt/alpha;
02771 }
02772
02773
02774
02775
02776
02777
02778
02779
02780 solution_correction ();
02781
02782 approximation ();
02783
02784
02785
02786 if (Tp->nvs==1 && Tp->pnvs==1)
02787 actual_previous_nodval ();
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797
02798
02799 }
02800
02801
02802 ani=j;
02803
02804 if (ani==0)
02805 nsts++;
02806
02807 if (ani==ini){
02808
02809 for (j=0;j<n;j++){
02810 lhs[j]=lhsb[j];
02811 tdlhs[j]=tdlhsb[j];
02812 }
02813
02814
02815
02816
02817
02818 solution_correction ();
02819
02820 approximation ();
02821
02822
02823
02824
02825 if (Tp->nvs==1 && Tp->pnvs==1)
02826 actual_previous_nodval ();
02827
02828
02829
02830
02831
02832 dt/=2.0;
02833 i--;
02834 Tp->timecont.oldtime ();
02835
02836
02837 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium dt = %le",dt);
02838
02839 if (dt<dtmin){
02840 fprintf (stderr,"\n\n time increment is less than minimum time increment");
02841 fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
02842 abort ();
02843 }
02844 }else{
02845
02846
02847
02848
02849 if (nsts==2){
02850 dt*=2.0;
02851
02852 nsts=0;
02853
02854
02855 fprintf (stdout,"\n time increment is enlarged because no inner loop was neccessary in previous 3 steps dt = %le",dt);
02856
02857 }
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870 print_stept(0,i,Tp->time,NULL);
02871 print_flusht();
02872 }
02873
02874 }while(Tp->time<end_time);
02875
02876 delete [] v;
02877 delete [] f;
02878 delete [] p;
02879 delete [] d;
02880 delete [] lhsb;
02881 delete [] tdlhsb;
02882
02883 print_closet();
02884 }
02885
02886
02887
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899
02900 void nonlinear_nonstat_solv_nr_dform ()
02901 {
02902 long i,j,k,n,ini,ani,nsts,lcid;
02903 long stop;
02904 double dt,dtmin,end_time,time,alpha,*err,norrhs,zero,*thresh;
02905 double *d,*f,*p,*v,*lhs,*lhsi,*tdlhs,*rhs,*lhsb,*tdlhsb;
02906 double *z,*diag;
02907
02908
02909 n=Ndoft;
02910
02911 lcid=0;
02912
02913 zero=1.0e-20;
02914
02915 thresh=Tp->threshrhs;
02916
02917
02918 lhs = Lsrst->give_lhs (0);
02919
02920 tdlhs = Lsrst->give_tdlhs (0);
02921
02922 lhsi = Lsrst->lhsi;
02923
02924 rhs = Lsrst->give_rhs (0);
02925
02926
02927
02928 d = new double [n];
02929
02930 p = new double [n];
02931
02932 f = new double [n];
02933
02934 v = new double [n];
02935
02936 lhsb = new double [n];
02937
02938 tdlhsb = new double [n];
02939
02940 z = new double [n];
02941 diag = new double [n];
02942
02943
02944
02945 nullv (lhs,n);
02946 nullv (rhs,n);
02947 nullv (tdlhs,n);
02948 nullv (d,n);
02949 nullv (p,n);
02950 nullv (f,n);
02951 nullv (v,n);
02952 nullv (z,n);
02953
02954
02955 approximation ();
02956
02957
02958 alpha=Tp->alpha;
02959
02960 err=Tp->errarr;
02961
02962 ini=Tp->nii;
02963
02964 dtmin=Tp->timecont.dtmin;
02965
02966
02967 Tp->time=Tp->timecont.starttime ();
02968 time=Tp->time;
02969
02970 end_time = Tp->timecont.endtime ();
02971
02972 dt = Tp->timecont.initialtimeincr ();
02973
02974
02975 Tm->initmaterialmodels();
02976
02977 approximation ();
02978
02979
02980
02981 if (Tp->nvs==1 && Tp->pnvs==1)
02982 actual_previous_nodval ();
02983
02984
02985
02986
02987
02988
02989
02990
02991 i=0;
02992
02993 nsts=0;
02994
02995 stop=0;
02996
02997
02998 Tm->initmaterialmodels();
02999 print_initt(-1, "wt");
03000 print_stept(0,i,Tp->time,NULL);
03001
03002
03003
03004 do{
03005
03006 Tp->time=Tp->timecont.newtime (dt);
03007
03008 if (Mesprt != 0) fprintf (stdout,"\n\n time %le, time step = %le, increment number %ld",Tp->time,dt,i);
03009
03010 Tm->updateipval ();
03011 i++;
03012
03013
03014 for (j=0;j<n;j++){
03015 lhsb[j]=lhs[j];
03016 tdlhsb[j]=tdlhs[j];
03017 }
03018
03019
03020 capacity_matrix (lcid);
03021
03022
03023
03024
03025
03026 conductivity_matrix (lcid);
03027
03028
03029
03030
03031
03032
03033
03034
03035
03036
03037 for (j=0;j<n;j++){
03038 d[j] = lhs[j]+dt*(1.0-alpha)*tdlhs[j];
03039 }
03040
03041
03042 Cmat->gmxv (d,p);
03043
03044
03045
03046 Kmat->scalgm (dt*alpha);
03047 Kmat->addgm (1.0,*Cmat);
03048
03049
03050
03051 trfel_right_hand_side (lcid,f,n);
03052
03053
03054 norrhs=0.0;
03055 for (j=0;j<n;j++){
03056 rhs[j] = f[j]*alpha*dt + p[j];
03057
03058
03059
03060
03061 }
03062
03063
03064
03065
03066
03067
03068
03069
03070
03071 Kmat->diag_check (zero,rhs);
03072
03073
03074
03075
03076
03077
03078
03079
03080
03081
03082 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
03083
03084
03085
03086
03087
03088
03089
03090
03091 for (j=0;j<n;j++){
03092 tdlhs[j]=(lhs[j]-d[j])/dt/alpha;
03093 }
03094
03095
03096
03097
03098
03099
03100
03101 solution_correction ();
03102
03103 approximation ();
03104
03105
03106
03107 if (Tp->nvs==1 && Tp->pnvs==1)
03108 actual_previous_nodval ();
03109
03110
03111
03112 for (j=0;j<ini;j++){
03113
03114
03115
03116
03117
03118
03119
03120
03121
03122
03123
03124
03125
03126
03127 Kmat->gmxv (lhs,v);
03128
03129
03130
03131
03132
03133
03134
03135 for (k=0;k<n;k++){
03136 rhs[k]=f[k]*alpha*dt+p[k]-v[k];
03137
03138
03139 z[k]=f[k]*alpha*dt+p[k];
03140
03141
03142 }
03143
03144
03145
03146
03147
03148 stop = norm_computation_vec (rhs,z,err,thresh,2,1);
03149
03150
03151
03152 if (stop==1){
03153 break;
03154 }
03155
03156
03157
03158
03159
03160
03161
03162
03163
03164
03165
03166
03167
03168
03169
03170
03171
03172
03173
03174
03175
03176
03177
03178
03179
03180
03181
03182
03183 if (Mesprt != 0) fprintf (stdout,"\n iteration number %ld",j);
03184
03185
03186 for (k=0;k<n;k++){
03187 rhs[k]=f[k]*alpha*dt+p[k];
03188
03189 }
03190
03191
03192
03193
03194
03195
03196 Kmat->diag_check (zero,rhs);
03197
03198
03199
03200
03201
03202
03203
03204
03205 Tp->ssle->solve_system (Gtt,Kmat,lhs,rhs,Outt);
03206
03207
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226
03227 for (k=0;k<n;k++){
03228 tdlhs[k]=(lhs[k]-d[k])/dt/alpha;
03229 }
03230
03231
03232
03233
03234
03235
03236
03237
03238 solution_correction ();
03239
03240 approximation ();
03241
03242
03243
03244 if (Tp->nvs==1 && Tp->pnvs==1)
03245 actual_previous_nodval ();
03246
03247
03248
03249
03250
03251
03252
03253
03254
03255
03256
03257 }
03258
03259
03260 ani=j;
03261
03262 if (ani==0)
03263 nsts++;
03264
03265 if (ani==ini){
03266
03267 for (j=0;j<n;j++){
03268 lhs[j]=lhsb[j];
03269 tdlhs[j]=tdlhsb[j];
03270 }
03271
03272
03273
03274
03275
03276 solution_correction ();
03277
03278 approximation ();
03279
03280
03281
03282
03283 if (Tp->nvs==1 && Tp->pnvs==1)
03284 actual_previous_nodval ();
03285
03286
03287
03288
03289
03290 dt/=2.0;
03291 i--;
03292 Tp->timecont.oldtime ();
03293
03294
03295 fprintf (stdout,"\n time increment is reduced because the inner loop was not able to enforce equilibrium dt = %le",dt);
03296
03297 if (dt<dtmin){
03298 fprintf (stderr,"\n\n time increment is less than minimum time increment");
03299 fprintf (stderr,"\n computation fails (file %s, line %d)\n",__FILE__,__LINE__);
03300 abort ();
03301 }
03302 }else{
03303
03304
03305
03306
03307 if (nsts==2){
03308 dt*=2.0;
03309
03310 nsts=0;
03311
03312
03313 fprintf (stdout,"\n time increment is enlarged because no inner loop was neccessary in previous 3 steps dt = %le",dt);
03314
03315 }
03316
03317
03318
03319
03320
03321
03322
03323
03324
03325
03326
03327
03328 print_stept(0,i,Tp->time,NULL);
03329 print_flusht();
03330 }
03331
03332 }while(Tp->time<end_time);
03333
03334 delete [] v;
03335 delete [] f;
03336 delete [] p;
03337 delete [] d;
03338 delete [] lhsb;
03339 delete [] tdlhsb;
03340
03341 print_closet();
03342 }