00001 #include <math.h>
00002 #include <string.h>
00003 #include "fcsolver.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 "gmatrix.h"
00011 #include "gtopology.h"
00012 #include "elemswitch.h"
00013 #include "dloadcase.h"
00014 #include "transprint.h"
00015 #include "mechprint.h"
00016 #include "backupsol.h"
00017 #include "backupsolt.h"
00018
00019
00020
00021
00022
00023
00024
00025 void solve_fcouplprob ()
00026 {
00027
00028
00029 double *rhs;
00030
00031 Lsrs->lhs=Lsrsc->lhs;
00032 Lsrs->rhs=Lsrsc->rhs;
00033 Lsrst->lhs=Lsrsc->lhs;
00034 Lsrst->rhs=Lsrsc->rhs;
00035
00036 rhs = Lsrsc->rhs;
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 nonlinear_solver_coupl (0);
00047
00048 Lsrs->lhs=NULL;
00049 Lsrs->rhs=NULL;
00050 Lsrst->lhs=NULL;
00051 Lsrst->rhs=NULL;
00052 }
00053
00054
00055 void nonlinear_solver_coupl (long lcid)
00056 {
00057 switch (Cp->tnlinsol){
00058 case newtonc:{
00059
00060 newton_raphson_coupl_new (lcid);
00061 break;
00062 }
00063 default:{
00064 fprintf (stderr,"\n\n unknown solver of nonlinear equation system is required");
00065 fprintf (stderr,"\n in function nonlinear_solver (file %s, line %d).\n",__FILE__,__LINE__);
00066 }
00067 }
00068
00069 }
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 void newton_raphson_coupl_new (long lcid)
00081 {
00082 long i,j,k,n,ini;
00083 double zero,dt,end_time,alpha,norfb,err;
00084 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*arhs;
00085 double *fp,*fl;
00086
00087
00088 alpha=Cp->alpha;
00089 zero=Cp->zero;
00090
00091 ini=Cp->niilnr;
00092
00093 err=Cp->errnr;
00094
00095
00096
00097 n=Ndofc;
00098
00099
00100
00101
00102
00103 lhs = Lsrsc->give_lhs (lcid);
00104
00105 tdlhs = Lsrsc->give_tdlhs (lcid);
00106
00107 lhsi = Lsrsc->lhsi;
00108
00109 rhs = Lsrsc->give_rhs (lcid);
00110
00111
00112 nullv (lhs,n);
00113 nullv (tdlhs,n);
00114 nullv (rhs,n);
00115
00116
00117
00118
00119
00120 Lsrs->lhs = lhs;
00121 Lsrs->rhs = rhs;
00122
00123
00124
00125 Lsrst->lhs = lhs;
00126 Lsrst->tdlhs = tdlhs;
00127 Lsrst->rhs = rhs;
00128
00129
00130
00131 d = new double [n];
00132 nullv (d,n);
00133
00134 p = new double [n];
00135 nullv (p,n);
00136
00137 arhs = new double [n];
00138 nullv (arhs,n);
00139
00140 fp = new double [n];
00141 nullv (fp,n);
00142
00143 fb = new double [n];
00144 nullv (fb,n);
00145
00146 fi = new double [n];
00147 nullv (fi,n);
00148 fl = new double [n];
00149 nullv (fi,n);
00150
00151
00152
00153 Cp->time=Cp->timecon.starttime ();
00154
00155 dt=Cp->timecon.initialtimeincr ();
00156
00157 end_time = Cp->timecon.endtime ();
00158
00159 Mp->time = Cp->time;
00160 Mp->timecon.take_values (Cp->timecon);
00161
00162 Tp->time = Cp->time;
00163 Tp->timecont.take_values (Cp->timecon);
00164
00165
00166
00167 approximationc ();
00168
00169
00170 Cmu->initmaterialmodels();
00171
00172
00173
00174
00175
00176 i=0;
00177
00178 if (Mp->hdbcont.restore_stat()){
00179 if (Mespr==1)
00180 fprintf (stdout,"\n Reading of backup file for MEFEL\n");
00181
00182
00183
00184 init_trfel_mefel ();
00185 Mm->initmaterialmodels(lcid);
00186 print_init(-1, "at");
00187 }
00188 else
00189 {
00190
00191
00192 init_trfel_mefel ();
00193 Mm->initmaterialmodels(lcid);
00194
00195 if (Mp->eigstrains > 0)
00196 internal_forces(lcid, fp);
00197 print_init(-1, "wt");
00198 print_step(lcid, i, Cp->time, rhs);
00199 print_flush();
00200 }
00201 if (Tp->hdbcont.restore_stat()){
00202 if (Mesprt==1)
00203 fprintf (stdout,"\n Reading of backup file for TRFEL\n");
00204
00205
00206 Tm->initmaterialmodels();
00207 compute_req_valt (0);
00208 print_initt(-1, "at");
00209 }
00210 else{
00211
00212 Tm->initmaterialmodels();
00213 compute_req_valt (0);
00214 print_initt(-1, "wt");
00215 print_stept(0,i,Tp->time,NULL);
00216 print_flusht();
00217 }
00218
00219 do{
00220
00221
00222 Cp->time = Cp->timecon.newtime ();
00223 dt = Cp->timecon.actualbacktimeincr ();
00224 Mp->time = Cp->time;
00225 Mp->timecon.take_values (Cp->timecon);
00226 Tp->time = Cp->time;
00227 Tp->timecont.take_values (Cp->timecon);
00228
00229 if (Mesprc==1) fprintf (stdout,"\n\n -------------------------------------------------------------------------------------------");
00230 if (Mesprc==1) fprintf (stdout,"\n NEW METR TIME STEP = %ld: METR TIME = %e, metr time increment = %e",i,Cp->time,dt);
00231 if (Mesprc==1) fprintf (stdout,"\n -------------------------------------------------------------------------------------------");
00232
00233
00234 approximationc ();
00235
00236 updateval();
00237
00238 copy_data();
00239
00240 i++;
00241
00242
00243 zero_order_matrix (lcid);
00244
00245
00246 first_order_matrix (lcid);
00247
00248
00249 for (j=0;j<n;j++){
00250 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00251 }
00252
00253
00254 metr_right_hand_side (lcid,rhs,fl);
00255
00256
00257 D0mat->gmxv (d,p);
00258
00259
00260
00261 D0mat->scalgm (dt*alpha);
00262 D0mat->addgm (1.0,*D1mat);
00263 D0mat->copygm (*D1mat);
00264
00265
00266
00267
00268
00269
00270 for (j=0;j<n;j++){
00271 fb[j] = rhs[j]-p[j];
00272 }
00273
00274
00275 Cp->ssle->solve_system (Gtu,D1mat,tdlhs,fb,Outc);
00276
00277
00278
00279 for (j=0;j<n;j++){
00280 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00281 }
00282
00283
00284 for (j=0;j<ini;j++){
00285
00286 solution_correction ();
00287
00288 approximationc ();
00289
00290 updateval();
00291
00292 copy_data();
00293
00294 if (Cp->fcsolv == fullnewtonc){
00295
00296
00297 if (D0mat != NULL){
00298 delete D0mat;
00299 D0mat=NULL;
00300 }
00301
00302 if (D1mat != NULL){
00303 delete D1mat;
00304 D1mat=NULL;
00305 }
00306
00307
00308 zero_order_matrix (lcid);
00309
00310
00311 first_order_matrix (lcid);
00312
00313
00314 D0mat->gmxv (d,p);
00315
00316
00317
00318 D0mat->scalgm (dt*alpha);
00319 D0mat->addgm (1.0,*D1mat);
00320 D0mat->copygm (*D1mat);
00321 }
00322
00323
00324 if (Cp->restype == fluxesc){
00325
00326
00327 internal_gforces (lcid,fi);
00328 for (k=0;k<n;k++){
00329 fb[k] = fl[k] - fi[k];
00330 }
00331 }
00332 if (Cp->restype == lrhsc){
00333
00334 metr_right_hand_side (lcid,rhs,fl);
00335
00336 D0mat->gmxv (tdlhs,fp);
00337
00338 internal_gforces (lcid,fi);
00339
00340 for (k=0;k<n;k++){
00341 fb[k] = rhs[k] - p[k] - fp[k];
00342 }
00343 }
00344
00345 norfb = ss (fb,fb,n);
00346
00347
00348 fprintf (stdout,"\n Inner loop j = %ld norf = %e\n",j,norfb);
00349
00350 if (norfb<err) break;
00351
00352 Cp->ssle->solve_system (Gtu,D1mat,p,fb,Outc);
00353
00354 for (k=0;k<n;k++){
00355 tdlhs[k]+=p[k];
00356 lhs[k]+=alpha*dt*p[k];
00357 }
00358 }
00359
00360 compute_req_val (0);
00361 print_step(0,i,Cp->time,NULL);
00362 print_flush();
00363 compute_req_valt (0);
00364 print_stept(0,i,Cp->time,NULL);
00365 print_flusht();
00366
00367 }while(Cp->time<end_time);
00368
00369 print_close ();
00370 print_closet ();
00371
00372 delete [] fp;
00373 delete [] fl;
00374 delete [] fi;
00375 delete [] fb;
00376 delete [] p;
00377 delete [] arhs;
00378 delete [] d;
00379 }
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390 void newton_raphson_coupl (long lcid)
00391 {
00392 long i,j,k,n,ini;
00393 double zero,dt,end_time,alpha,norfb,err;
00394 double *d,*fb,*fi,*p,*lhs,*lhsi,*tdlhs,*rhs,*arhs;
00395
00396
00397 alpha=Cp->alpha;
00398 zero=Cp->zero;
00399
00400 ini=Cp->niilnr;
00401
00402 err=Cp->errnr;
00403
00404
00405
00406 n=Ndofc;
00407
00408
00409
00410
00411
00412 lhs = Lsrsc->give_lhs (lcid);
00413
00414 tdlhs = Lsrsc->give_tdlhs (lcid);
00415
00416 lhsi = Lsrsc->lhsi;
00417
00418 rhs = Lsrsc->give_rhs (lcid);
00419
00420
00421 nullv (lhs,n);
00422 nullv (tdlhs,n);
00423 nullv (rhs,n);
00424
00425
00426
00427
00428
00429 Lsrs->lhs = lhs;
00430 Lsrs->rhs = rhs;
00431
00432
00433
00434 Lsrst->lhs = lhs;
00435 Lsrst->tdlhs = tdlhs;
00436 Lsrst->rhs = rhs;
00437
00438
00439
00440 d = new double [n];
00441 nullv (d,n);
00442
00443 p = new double [n];
00444 nullv (p,n);
00445
00446 arhs = new double [n];
00447 nullv (arhs,n);
00448
00449 fb = new double [n];
00450 nullv (fb,n);
00451
00452 fi = new double [n];
00453 nullv (fi,n);
00454
00455
00456 Cp->time=Cp->timecon.starttime ();
00457
00458 dt=Cp->timecon.initialtimeincr ();
00459
00460 end_time = Cp->timecon.endtime ();
00461
00462
00463 approximationc ();
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473 compute_req_valt (0);
00474
00475
00476
00477
00478 i=0;
00479 print_init(-1, "wt");
00480 print_step(0,i,Cp->time,NULL);
00481 print_initt(-1, "wt");
00482 print_stept(0,i,Cp->time,NULL);
00483 do{
00484
00485 if (Mesprc==1) fprintf (stdout,"\n time %e",Cp->time);
00486
00487
00488 zero_order_matrix (lcid);
00489
00490
00491 first_order_matrix (lcid);
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 for (j=0;j<n;j++){
00503 d[j] = lhs[j] + (1.0-alpha)*dt*tdlhs[j];
00504 }
00505
00506
00507 metr_right_hand_side (lcid,rhs,p);
00508
00509
00510 D0mat->gmxv (d,p);
00511
00512
00513
00514 D0mat->scalgm (dt*alpha);
00515 D0mat->addgm (1.0,*D1mat);
00516
00517
00518 for (j=0;j<n;j++){
00519 fb[j] = rhs[j] - p[j];
00520 }
00521
00522
00523
00524
00525
00526
00527
00528
00529 Cp->ssle->solve_system (Gtt,D0mat,tdlhs,fb,Outt);
00530
00531
00532
00533
00534
00535
00536
00537 for (j=0;j<n;j++){
00538 lhs[j] = d[j] + alpha*dt*tdlhs[j];
00539 }
00540
00541
00542
00543
00544
00545
00546 solution_correction ();
00547 approximationc ();
00548
00549
00550
00551 internal_gforces (lcid,fi);
00552
00553
00554
00555
00556
00557
00558
00559 if (Cp->restype == fluxesc){
00560 for (k=0;k<n;k++){
00561 fb[k]=fi[k];
00562 }
00563 }
00564 if (Cp->restype == lrhsc){
00565 D0mat->gmxv (tdlhs,fi);
00566 for (k=0;k<n;k++){
00567 fb[k] = rhs[k] - p[k] - fi[k];
00568 }
00569 }
00570
00571 norfb = ss (fb,fb,n);
00572
00573 if (Mesprc==1) fprintf (stdout,"\n%e %e",norfb,err);
00574
00575 if (norfb<err){
00576 Cp->time = Cp->timecon.newtime ();
00577 dt = Cp->timecon.actualbacktimeincr ();
00578 i++;
00579
00580 Mp->time = Cp->time;
00581 Mp->timecon.take_values (Cp->timecon);
00582
00583
00584 Tp->time = Cp->time;
00585 Tp->timecont.take_values (Cp->timecon);
00586
00587 print_step(0,i,Cp->time,NULL);
00588 print_flush();
00589 print_stept(0,i,Cp->time,NULL);
00590 print_flusht();
00591 continue;
00592 }
00593
00594
00595 for (j=0;j<ini;j++){
00596
00597 if (Cp->fcsolv == fullnewtonc){
00598
00599 zero_order_matrix (lcid);
00600
00601
00602 first_order_matrix (lcid);
00603
00604
00605 D0mat->gmxv (d,p);
00606
00607
00608
00609 D0mat->scalgm (dt*alpha);
00610 D0mat->addgm (1.0,*D1mat);
00611 }
00612
00613
00614 Cp->ssle->solve_system (Gtt,D0mat,p,fb,Outt);
00615
00616 for (k=0;k<n;k++){
00617 tdlhs[k]+=p[k];
00618 lhs[k]+=alpha*dt*p[k];
00619 }
00620
00621 solution_correction ();
00622 approximationc ();
00623 internal_gforces (lcid,fi);
00624
00625
00626
00627
00628
00629
00630
00631 if (Cp->restype == fluxesc){
00632 for (k=0;k<n;k++){
00633 fb[k]=fi[k];
00634 }
00635 }
00636 if (Cp->restype == lrhsc){
00637 D0mat->gmxv (tdlhs,fi);
00638 for (k=0;k<n;k++){
00639 fb[k] = rhs[k] - p[k] - fi[k];
00640 }
00641 }
00642
00643 norfb = ss (fb,fb,n);
00644
00645 if (Mesprc==1) fprintf (stdout,"\n iterace %ld chyba %e",j,norfb);
00646
00647 if (norfb<err){
00648 break;
00649 }
00650 }
00651
00652
00653
00654
00655 Cp->time = Cp->timecon.newtime ();
00656 dt = Cp->timecon.actualbacktimeincr ();
00657 i++;
00658
00659 Mp->time = Cp->time;
00660 Mp->timecon.take_values (Cp->timecon);
00661
00662
00663 Tp->time = Cp->time;
00664 Tp->timecont.take_values (Cp->timecon);
00665
00666 print_step(0,i,Cp->time,NULL);
00667 print_flush();
00668 print_stept(0,i,Cp->time,NULL);
00669 print_flusht();
00670
00671 fprintf (Outc,"\n\n\n posun konce tyce d %f",lhs[Gtm->gnodes[4].cn[0]]);
00672
00673 }while(Cp->time<end_time);
00674
00675
00676 delete [] fi;
00677 delete [] fb;
00678 delete [] p;
00679 delete [] arhs;
00680 delete [] d;
00681 }
00682
00683
00684 void vector_assemb (double *c,double *m,double *t)
00685 {
00686 long i,j;
00687
00688 for (i=0;i<Gtm->ndof;i++){
00689 j=Gtm->cngtopcorr[i];
00690 c[j]=m[i];
00691 }
00692 for (i=0;i<Gtt->ndof;i++){
00693 j=Gtt->cngtopcorr[i];
00694 c[j]=t[i];
00695 }
00696
00697 }
00698
00699 void vector_decomp (double *c,double *m,double *t)
00700 {
00701 long i,j;
00702
00703 for (i=0;i<Gtm->ndof;i++){
00704 j=Gtm->cngtopcorr[i];
00705 m[i]=c[j];
00706 }
00707 for (i=0;i<Gtt->ndof;i++){
00708 j=Gtt->cngtopcorr[i];
00709 t[i]=c[j];
00710 }
00711
00712 }