00001 #include <string.h>
00002 #include "mtsolver.h"
00003 #include "backupsol.h"
00004 #include "mtglvec.h"
00005 #include "global.h"
00006 #include "globmat.h"
00007 #include "mechprint.h"
00008 #include "matrix.h"
00009 #include "vector.h"
00010 #include "gmatrix.h"
00011 #include "gtopology.h"
00012 #include "dloadcase.h"
00013 #include "elemswitch.h"
00014 #include "nssolver.h"
00015 #include "newtonraph.h"
00016 #include "mathem.h"
00017
00018
00019 #ifndef FNAMELEN
00020 #define FNAMELEN 1001
00021 #endif
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 void solve_time_dep_prob ()
00034 {
00035 long lcid;
00036
00037
00038
00039 lcid=0;
00040
00041
00042 visco_solver2 (lcid);
00043 }
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void visco_solver (long lcid)
00058 {
00059 long i,j,k,n,ini,nsts;
00060 double dt,dtmin,dtdef,end_time;
00061 double *f,*fl,*fi,*fb,*fp,*r,*dr,*lhsb;
00062 double norf, norfa, err, zero;
00063 matrix lsm_a(3,3);
00064 vector lsm_r(3), lsm_l(3);
00065 Mp->nlman->divc_step = new double[10];
00066 Mp->nlman->divc_err = new double[10];
00067
00068
00069
00070
00071 ini = Mp->nlman->niilnr;
00072
00073 err = Mp->nlman->errnr;
00074
00075 n=Ndofm;
00076
00077 zero = Mp->zero;
00078
00079
00080 r = Lsrs->give_lhs (lcid);
00081 nullv (r,n);
00082
00083 f = Lsrs->give_rhs (lcid);
00084
00085
00086
00087
00088 dr = new double [n];
00089 nullv (dr,n);
00090
00091 fp = new double [n];
00092 nullv (fp,n);
00093
00094 fl = new double [n];
00095 nullv (fl,n);
00096
00097 fi = new double [n];
00098 nullv (fi,n);
00099
00100 fb = new double [n];
00101 nullv (fb,n);
00102
00103 lhsb = new double [n];
00104 nullv (lhsb,n);
00105
00106
00107
00108 Mp->time=Mp->timecon.starttime ();
00109
00110 dt = Mp->timecon.initialtimeincr ();
00111 Mp->dtime = dt;
00112
00113 end_time = Mp->timecon.endtime ();
00114
00115 dtmin=Mp->timecon.dtmin;
00116
00117
00118 i=0;
00119
00120 nsts=0;
00121
00122
00123 if (Mp->hdbcont.restore_stat()){
00124 if (Mespr==1)
00125 fprintf (stdout,"\n Reading of backup file\n");
00126 solver_restore (r,fp,i,Mp->time,dt,&Mp->timecon,n);
00127
00128 Mm->initmaterialmodels(lcid);
00129 print_init(-1, "at");
00130
00131
00132 }
00133 else
00134 {
00135
00136 Mm->initmaterialmodels(lcid);
00137 print_init(-1, "wt");
00138 print_step(lcid, i, Mp->time, f);
00139
00140 print_flush();
00141 }
00142
00143
00144
00145
00146
00147 do{
00148
00149 Mp->time = Mp->timecon.newtime (dt);
00150 Mp->dtime = dt;
00151
00152 i++;
00153 Mp->istep = i;
00154 Mp->jstep = -1;
00155
00156
00157 Mp->strainstate=0;
00158
00159
00160 Mp->stressstate=0;
00161
00162
00163 Mp->otherstate=0;
00164
00165
00166 for (j=0;j<n;j++){
00167 lhsb[j]=r[j];
00168 }
00169
00170 if (Mespr==1) fprintf (stdout,"\n\n -------------------------------------------------------------");
00171 if (Mespr==1) fprintf (stdout,"\n MEFEL Time step = %ld, Time %e, Time increment = %e",i,Mp->time,dt);
00172 if (Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------\n");
00173
00174
00175 mefel_right_hand_side (lcid,f,fl);
00176
00177
00178
00179 Mb->comp_sum (f);
00180
00181
00182
00183
00184
00185 incr_internal_forces (lcid,fi);
00186
00187
00188 for (j=0;j<n;j++){
00189 fb[j]=f[j]-fp[j]+fi[j];
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 if ((Mp->nlman->stmat==tangent_stiff) || (i == 1) || (Smat==NULL))
00203 stiffness_matrix (lcid);
00204
00205
00206 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00207
00208
00209 for (j=0;j<n;j++){
00210 r[j]+=dr[j];
00211 }
00212
00213
00214 internal_forces (lcid,fi);
00215
00216
00217 norfa=ss(f,f,n);
00218
00219 for (j=0;j<n;j++){
00220 fb[j] = fl[j] - fi[j];
00221 }
00222
00223
00224 norf=ss(fb,fb,n);
00225 if (norfa > 0.0)
00226 norf /= norfa;
00227
00228 if (Mespr==1) fprintf (stdout,"\n\n Norf before inner iteration loop norf = %e\n\n\n",norf);
00229 j = 0;
00230 if (norf<err){
00231 nsts++;
00232 }
00233 else{
00234
00235
00236 fillm(0.0, lsm_a);
00237 fillv(0.0, lsm_r);
00238 for (j=0;j<ini;j++)
00239 {
00240 Mp->jstep = j;
00241
00242 if ((Mp->nlman->stmat==tangent_stiff) && (j%5 == 0))
00243 stiffness_matrix (lcid);
00244
00245 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00246
00247
00248 for (k=0;k<n;k++)
00249 r[k]+=dr[k];
00250
00251
00252 internal_forces (lcid,fi);
00253
00254
00255 for (k=0;k<n;k++)
00256 fb[k]=fl[k] - fi[k];
00257
00258
00259 norf=ss(fb,fb,n);
00260 if (norfa > 0.0)
00261 norf /= norfa;
00262
00263
00264 if (Mespr==1)
00265 fprintf (stdout,"\n Inner loop number j = %ld, norf = %e \n",j,norf);
00266
00267 if (norf<err) break;
00268
00269
00270 if (j > 10)
00271 {
00272 if (lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), Mp->nlman->divc_step[j%10],
00273 Mp->nlman->divc_err[j%10], norf, Mp->zero,1) > 0.0)
00274 {
00275 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00276 break;
00277 }
00278 Mp->nlman->divc_step[j%10] = double(j+1);
00279 Mp->nlman->divc_err[j%10] = err;
00280 }
00281 else
00282 {
00283 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
00284 Mp->nlman->divc_step[j%10] = double(j+1);
00285 Mp->nlman->divc_err[j%10] = err;
00286 }
00287 }
00288 }
00289
00290 if (j==ini || norf>err)
00291 {
00292
00293 if (Mespr==1) fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00294
00295
00296 for (j=0;j<n;j++)
00297 r[j]=lhsb[j];
00298
00299
00300
00301 dt/=2.0;
00302 Mp->timecon.oldtime ();
00303 i--;
00304
00305 if (Mespr==1) fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00306 if (dt<dtmin)
00307 {
00308 if (Mespr==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00309 if (Mespr==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00310 if (Mespr==1) fprintf (stderr,"\n FORCED output of results from this step is performed\n");
00311
00312 Mm->updateipval();
00313 compute_req_val (lcid);
00314 print_step_forced(lcid, i, Mp->time, fb);
00315
00316 print_flush();
00317 break;
00318 }
00319 }
00320 else
00321 {
00322
00323 if (Mespr==1) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00324
00325
00326 for (j=0;j<n;j++)
00327 fp[j]=f[j];
00328
00329 dtdef = Mp->timecon.actualforwtimeincr ();
00330 if (nsts==2)
00331 {
00332 dt*=2.0;
00333 nsts=0;
00334 if (Mespr==1) fprintf (stdout,"\n\n time increment is enlarged to %le because no inner loop was neccessary in previous 3 steps", dt);
00335 }
00336
00337
00338
00339
00340 Mb->comp_sum (fi);
00341
00342
00343
00344
00345
00346 Mm->updateipval();
00347 compute_req_val (lcid);
00348 print_step(lcid, i, Mp->time, fi);
00349
00350 print_flush();
00351
00352 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.save_stat()))
00353 {
00354 if (Mespr==1)
00355 fprintf (stdout,"\n Creating MEFEL backup file\n");
00356 solver_save (r,fp,i,Mp->time,dt,&Mp->timecon,n);
00357 }
00358 }
00359 }while(Mp->time<=end_time);
00360
00361 print_close ();
00362
00363 delete [] fi;
00364 delete [] fb;
00365 delete [] fp;
00366 delete [] dr;
00367 delete [] fl;
00368 delete [] lhsb;
00369 }
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 void visco_solver2 (long lcid)
00384 {
00385 long i, li, nsts, ret;
00386 double time, dt, dtmin, dtmax, dtr, end_time;
00387 mt_glob_vec mt_gv;
00388
00389
00390
00391
00392 visco_solver_init(lcid, mt_gv);
00393
00394
00395 dt = Mp->timecon.initialtimeincr ();
00396
00397 dtmin=Mp->timecon.dtmin;
00398
00399 dtmax=Mp->timecon.dtmax;
00400
00401 end_time = Mp->timecon.endtime ();
00402
00403
00404 li = i = mt_gv.istep;
00405
00406 nsts=0;
00407
00408 ret = 0;
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424 do
00425 {
00426
00427 i++;
00428
00429
00430 time=Mp->timecon.newtime(dt);
00431 Mp->dtime = dt;
00432 ret = one_step(lcid, time, dt, dtr, i, li, mt_gv);
00433
00434 if (ret >= 0)
00435 {
00436
00437
00438 copyv(mt_gv.f, mt_gv.fp, Ndofm);
00439
00440
00441 Mm->updateipval();
00442 compute_req_val (lcid);
00443 print_step(lcid, Mp->istep, Mp->time, mt_gv.fl);
00444 print_flush();
00445
00446 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.save_stat()))
00447 {
00448 if (Mespr==1)
00449 fprintf (stdout,"\n Creating MEFEL backup file\n");
00450 solver_save (mt_gv.r,mt_gv.fp,Mp->istep,Mp->time,dt,&Mp->timecon,Ndofm);
00451 }
00452 if (ret == 0)
00453 nsts++;
00454 if (nsts==2)
00455 {
00456 dt*=2.0;
00457 nsts=0;
00458 if (Mespr==1)
00459 fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00460 }
00461
00462 if (dt>dtmax)
00463 dt = dtmax;
00464 }
00465 else
00466 {
00467
00468
00469 if (dtr < 1.0)
00470 dt *= dtr;
00471 else
00472 dt/=2.0;
00473 Mp->timecon.oldtime ();
00474 i--;
00475
00476 if (Mespr==1)
00477 fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00478 if (dt<dtmin)
00479 {
00480 if (Mespr==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00481 if (Mespr==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00482 if (Mespr==1) fprintf (stderr,"\n FORCED output of results from this step is performed\n");
00483
00484 Mm->updateipval();
00485 compute_req_val (lcid);
00486 print_step_forced(lcid, i, Mp->time, mt_gv.fb);
00487
00488 print_flush();
00489 break;
00490 }
00491 }
00492 }while(Mp->time<end_time);
00493
00494 print_close ();
00495 }
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505 void visco_solver_init(long lcid, mt_glob_vec &mt_gv)
00506 {
00507 long n;
00508 double dt;
00509
00510 n=Ndofm;
00511
00512 mt_gv.alloc(n);
00513
00514 mt_gv.r = Lsrs->give_lhs (lcid);
00515
00516 mt_gv.f = Lsrs->give_rhs (lcid);
00517
00518
00519 Mp->time=Mp->timecon.starttime ();
00520
00521
00522 if (Mp->hdbcont.restore_stat())
00523 {
00524 if (Mespr==1)
00525 fprintf (stdout,"\n Reading of backup file\n");
00526 solver_restore (mt_gv.r, mt_gv.fp, mt_gv.istep, Mp->time, dt, &Mp->timecon, n);
00527
00528 Mm->initmaterialmodels(lcid);
00529 print_init(-1, "at");
00530
00531
00532 }
00533 else
00534 {
00535 mt_gv.istep=0;
00536
00537 Mm->initmaterialmodels(lcid);
00538
00539 if (Mp->eigstrains > 0)
00540 internal_forces(lcid, mt_gv.fp);
00541
00542 print_init(-1, "wt");
00543 print_step(lcid, mt_gv.istep, Mp->time, mt_gv.f);
00544
00545 print_flush();
00546 }
00547 }
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 long one_step (long lcid,double time, double dt, double &dtr, long istep, long li, mt_glob_vec &mt_gv)
00571 {
00572 long j,n,ini;
00573 double *f,*fl,*fi,*fb,*fp,*r,*dr,*lhsb;
00574 double norf, err;
00575 matrix lsm_a(3,3);
00576 vector lsm_r(3), lsm_l(3);
00577 Mp->nlman->divc_step = new double[10];
00578 Mp->nlman->divc_err = new double[10];
00579
00580
00581 Mp->time = time;
00582
00583 Mp->istep = istep;
00584 Mp->jstep = -1;
00585
00586
00587 Mp->strainstate=0;
00588
00589
00590 Mp->stressstate=0;
00591
00592
00593 Mp->otherstate=0;
00594
00595
00596 r = mt_gv.r;
00597
00598 f = mt_gv.f;
00599
00600 dr = mt_gv.dr;
00601
00602 fp = mt_gv.fp;
00603
00604 fl = mt_gv.fl;
00605
00606 fi = mt_gv.fi;
00607
00608 fb = mt_gv.fb;
00609
00610 lhsb = mt_gv.lhsb;
00611
00612
00613 n=Ndofm;
00614
00615
00616 ini = Mp->nlman->niilnr;
00617
00618 err = Mp->nlman->errnr;
00619
00620
00621 for (j=0;j<n;j++)
00622 lhsb[j]=r[j];
00623
00624
00625 if (Mespr==1) fprintf (stdout,"\n\n ------------------------------------------------------------------------");
00626 if (Mespr==1) fprintf (stdout,"\n MEFEL Time step = %ld, Time %e, Time increment = %e",istep,Mp->time,dt);
00627 if (Mespr==1) fprintf (stdout,"\n ------------------------------------------------------------------------\n");
00628
00629 if (Mespr==1) fprintf (Out,"\n\n ------------------------------------------------------------------------");
00630 if (Mespr==1) fprintf (Out,"\n MEFEL Time step = %ld, Time %e, Time increment = %e",istep,Mp->time,dt);
00631 if (Mespr==1) fprintf (Out,"\n ------------------------------------------------------------------------\n");
00632
00633
00634 mefel_right_hand_side (lcid,f,fl);
00635
00636
00637
00638 Mb->comp_sum (f);
00639
00640
00641
00642
00643
00644 incr_internal_forces (lcid,fi);
00645
00646
00647
00648 subv(f, fp, fb, n);
00649 addv(fb, fi, fb, n);
00650
00651 j = 0;
00652 norf = gnewton_raphson_one_step(lcid, Mp->nlman, fl, r, fb, dr, fi, dtr, istep-1, j, li);
00653
00654 if (norf>err)
00655 {
00656
00657 if (Mespr==1) fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00658
00659
00660
00661 copyv(lhsb, r, n);
00662
00663 return -1;
00664 }
00665
00666
00667 if ((Mespr==1) && (j > 0)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00668
00669
00670 Mb->comp_sum (fi);
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681 if (j < 0)
00682 return 0;
00683
00684 return j+1;
00685 }