00001 #include <string.h>
00002 #include "cpsolver.h"
00003 #include "mtsolver.h"
00004 #include "nssolver.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 "node.h"
00014 #include "element.h"
00015 #include "slipsurf.h"
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 void solve_prob_constr_phases ()
00027 {
00028 long i=0;
00029 for (i=0; i<Lsrs->nlc; i++)
00030 solve_prob_constr_phases(i);
00031
00032 if (Mm->mcoul)
00033 {
00034 fprintf(stdout, "\n Mohr-Coulomb material detected -> starting detection of slip surfaces\n");
00035 long *plast_ip = new long[Mm->tnip];
00036 long nplz = detect_plastic_zones(plast_ip, 1.0e-3);
00037 fprintf(stdout, "\n There are detected %ld of plastic zones\n", nplz);
00038 for(i=0; i<nplz; i++)
00039 {
00040 long err;
00041 double a, b;
00042 err = approx_slip_surf(i+1, plast_ip, a, b);
00043 if (err)
00044 fprintf(stdout, " - Approximation of the %ld slip surface could not be calculated\n", i+1);
00045 else
00046 fprintf(stdout, " - Approximation of the %ld slip surface - a=%le, b=%le\n", i+1, a, b);
00047 }
00048 delete [] plast_ip;
00049 }
00050
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 void solve_prob_constr_phases (long lcid)
00067 {
00068 long i,j,k,n,ini,nsts,mncd,mnce,mnae,li;
00069 long nifn;
00070 long *ifn;
00071 double dt,dtmin,dtdef,end_time,prev_timem;
00072 double *f,*fl,*fi,*fb,*fp,*r,*dr,*lhsb;
00073 double norf, norfa, err, zero;
00074 matrix lsm_a(3,3);
00075 vector lsm_r(3), lsm_l(3);
00076
00077
00078 ifn = new long[Mt->nn];
00079 memset(ifn, 0, sizeof(*ifn)*Mt->nn);
00080 Mp->nlman->divc_step = new double[10];
00081 Mp->nlman->divc_err = new double[10];
00082
00083
00084
00085 ini = Mp->niilnr;
00086
00087 err = Mp->errnr;
00088
00089
00090
00091 n=Ndofm;
00092
00093 zero = Mp->zero;
00094
00095
00096
00097 norf = 0.0;
00098
00099
00100 r = Lsrs->give_lhs (lcid);
00101
00102 f = Lsrs->give_rhs (lcid);
00103
00104
00105
00106
00107 dr = new double [n];
00108 nullv (dr,n);
00109
00110 fp = new double [n];
00111 nullv (fp,n);
00112
00113 fl = new double [n];
00114 nullv (fl,n);
00115
00116 fi = new double [n];
00117 nullv (fi,n);
00118
00119 fb = new double [n];
00120 nullv (fb,n);
00121
00122 lhsb = new double [n];
00123 nullv (lhsb,n);
00124
00125
00126
00127 Mp->time=Mp->timecon.starttime ();
00128
00129 dt = Mp->timecon.initialtimeincr ();
00130
00131 end_time = Mp->timecon.endtime ();
00132
00133 dtmin=Mp->timecon.dtmin;
00134
00135 prev_timem = Mp->time;
00136
00137
00138 i=0;
00139
00140 nsts=0;
00141
00142
00143 if (Mp->hdbcont.restore_stat()){
00144 if (Mespr==1)
00145 fprintf (stdout,"\n Reading of backup file\n");
00146 solver_restore (r,fp,i,Mp->time,dt,&Mp->timecon,n);
00147 li = i;
00148
00149 prev_timem = Mp->time;
00150
00151
00152
00153 copyv(r, lhsb, n);
00154
00155
00156 Gtm->update_elem (Mp->time);
00157
00158 Gtm->update_auxinf();
00159
00160 Gtm->update_nodes ();
00161
00162 Gtm->update_dofs (Mp->time);
00163
00164 n = Ndofm = Gtm->codenum_generation(Out);
00165
00166 Mt->save_nodval(lcid);
00167
00168 Mt->comreac();
00169
00170 Mt->elemprescdisp();
00171
00172 Mm->initmaterialmodels(lcid);
00173
00174
00175
00176
00177 }
00178 else
00179 {
00180
00181
00182 li = 0;
00183 print_init(-1, "wt");
00184 print_step(lcid, i, Mp->time, f);
00185 print_flush();
00186 }
00187
00188
00189
00190 Mt->save_nodval(lcid);
00191
00192
00193
00194
00195 do{
00196
00197 Mp->time = Mp->timecon.newtime (dt);
00198 if (prev_timem == Mp->timecon.starttime())
00199 prev_timem = Mp->time;
00200
00201 i++;
00202 Mp->istep = i;
00203
00204
00205 Mp->strainstate=0;
00206
00207
00208 Mp->stressstate=0;
00209
00210
00211 Mp->otherstate=0;
00212
00213 if (Mespr==1) fprintf (stdout,"\n\n -------------------------------------------------------------");
00214 if (Mespr==1) fprintf (stdout,"\n Time step = %ld, Time %e, Time increment = %e",i,Mp->time,dt);
00215 if (Mespr==1) fprintf (stdout,"\n -------------------------------------------------------------\n");
00216
00217
00218
00219 mnce = Gtm->search_changed_elem(Mp->time, mnae);
00220
00221
00222 mncd = Gtm->search_changed_dofs(Mp->time,prev_timem);
00223
00224 nifn = 0;
00225 if (mnce!=0){
00226
00227
00228 nifn = Gtm->search_iface_nodes(ifn);
00229 }
00230
00231 if ((mnce!=0) || (mncd!=0))
00232 {
00233
00234
00235 Mt->save_node_inidispl(lcid, Mp->time, prev_timem);
00236
00237
00238 Gtm->switch_new_elem();
00239
00240
00241 Gtm->update_nodes();
00242
00243
00244
00245 Gtm->update_active_dofs(Mp->time);
00246
00247
00248 Ndofm = Gtm->codenum_generation(Out);
00249
00250
00251
00252 Mt->restore_nodval(r, Ndofm);
00253
00254
00255 Mm->initmaterialmodels(lcid);
00256
00257
00258
00259
00260
00261
00262 if (nifn)
00263 {
00264 if (Mp->comp_inidispl)
00265 {
00266
00267 stress_initdispl(lcid);
00268
00269
00270
00271 Gtm->update_active_dofs(Mp->time);
00272
00273 Gtm->clear_intf_dofs(ifn);
00274
00275
00276 Ndofm = Gtm->codenum_generation(Out);
00277
00278
00279 Mp->strcomp = 0;
00280 internal_forces(lcid,fb);
00281 cmulv(-1.0, fb, Ndofm);
00282 Mp->strcomp = 1;
00283
00284
00285 delete Smat;
00286 Smat=NULL;
00287 stiffness_matrix (lcid);
00288 nullv(r, Ndofm);
00289
00290
00291
00292 Mp->ssle->solve_system (Gtm,Smat,r,fb,Out);
00293
00294
00295
00296
00297
00298
00299 Mt->save_elem_inidispl(lcid, ifn);
00300
00301
00302
00303 Gtm->remove_nodes(ifn);
00304
00305
00306
00307 Mt->save_nodval(lcid);
00308 }
00309 else
00310 {
00311
00312
00313
00314
00315
00316
00317 Mt->save_elem_inidispl(lcid, ifn);
00318 }
00319 }
00320
00321
00322 Gtm->update_elem(Mp->time);
00323
00324
00325 Gtm->update_nodes();
00326
00327
00328 Gtm->update_dofs (Mp->time);
00329
00330
00331 n = Ndofm = Gtm->codenum_generation(Out);
00332
00333
00334 delete Smat;
00335 Smat=NULL;
00336
00337
00338
00339 Mt->restore_nodval(r,n);
00340
00341
00342
00343
00344 if (i > 1)
00345 Mt->restore_nodforce(fp, n);
00346
00347
00348 Mt->comreac();
00349
00350 Mt->elemprescdisp();
00351
00352 Mt->clean_ip_new_elem();
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363 }
00364
00365
00366 mefel_right_hand_side (lcid,f,fl);
00367
00368
00369
00370
00371 Mb->comp_sum (f);
00372 Mb->comp_sum_react();
00373 Mb->comp_sum_pdreact();
00374 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], fx=%le, fy=%le, fz=%le",
00375 i, Mp->time, Mp->time/86400, Mb->sumcomp[0], Mb->sumcomp[1], Mb->sumcomp[2]);
00376 fflush(Out);
00377
00378 if ((mnce - mnae) > 0)
00379 {
00380
00381
00382
00383
00384 Mp->strcomp = 0;
00385 internal_forces(lcid,fi);
00386 Mp->strcomp = 1;
00387
00388
00389 subv(fp, fi, fb, n);
00390 }
00391 else
00392 nullv(fb, n);
00393
00394
00395 incr_internal_forces (lcid,fi);
00396
00397
00398 for (j=0;j<n;j++){
00399 fb[j]+=f[j]-fp[j]+fi[j];
00400 }
00401
00402
00403
00404
00405 assemble_stiffness_matrix(lcid, i, -1, li);
00406
00407
00408 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00409
00410
00411 addv(r, dr, n);
00412
00413
00414 internal_forces (lcid,fi);
00415
00416
00417 norfa=ss(f,f,n);
00418
00419
00420 subv(fl, fi, fb, n);
00421
00422
00423 norf=ss(fb,fb,n);
00424 if (norfa > 0.0)
00425 norf /= norfa;
00426
00427 if (Mespr==1) fprintf (stdout,"\n\n Norf before inner iteration loop norf = %e\n\n\n",norf);
00428 j = -1;
00429 if (norf<err){
00430 nsts++;
00431 }
00432 else
00433 {
00434
00435
00436 fillm(0.0, lsm_a);
00437 fillv(0.0, lsm_r);
00438 for (j=0;j<ini;j++)
00439 {
00440 Mp->jstep = j;
00441
00442
00443 assemble_stiffness_matrix(lcid,i,j+1,li);
00444
00445
00446 Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
00447
00448
00449 for (k=0;k<n;k++)
00450 r[k]+=dr[k];
00451
00452
00453 internal_forces (lcid,fi);
00454
00455
00456 subv(fl, fi, fb, n);
00457
00458
00459 norf=ss(fb,fb,n);
00460 if (norfa > 0.0)
00461 norf /= norfa;
00462
00463
00464 if (Mespr==1)
00465 fprintf (stdout,"\n Inner loop number j = %ld, norf = %e \n",j,norf);
00466
00467 if (norf<err) break;
00468
00469
00470 if (j > 10)
00471 {
00472 if (lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), Mp->nlman->divc_step[j%10],
00473 Mp->nlman->divc_err[j%10], norf, Mp->zero,1) > 0.0)
00474 {
00475 fprintf (stdout,"\n\n Divergence of inner iteation loop is detected. Inner loop is stopped.\n\n");
00476 break;
00477 }
00478 Mp->nlman->divc_step[j%10] = double(j+1);
00479 Mp->nlman->divc_err[j%10] = err;
00480 }
00481 else
00482 {
00483 lsm_quad(lsm_a, lsm_r, lsm_l, double(j+1), norf, 0.0, 0.0, Mp->zero,0);
00484 Mp->nlman->divc_step[j%10] = double(j+1);
00485 Mp->nlman->divc_err[j%10] = err;
00486 }
00487 }
00488 }
00489
00490 if (j==ini || norf>err)
00491 {
00492
00493 if (Mespr==1) fprintf (stdout,"\n\n Iteration process does not converge req. err = %e, reached err = %e", err,norf);
00494
00495
00496
00497
00498 dt/=2.0;
00499 Mp->timecon.oldtime ();
00500 i--;
00501
00502
00503
00504
00505 if ((mnce!=0) || (mncd!=0))
00506 {
00507
00508 Gtm->update_elem(prev_timem);
00509
00510 Gtm->update_nodes();
00511
00512 Gtm->update_dofs (prev_timem);
00513
00514 n = Ndofm = Gtm->codenum_generation(Out);
00515 }
00516
00517
00518 for (j=0;j<n;j++)
00519 r[j]=lhsb[j];
00520
00521 if (Mespr==1) fprintf (stdout,"\n\n time increment is reduced because the inner loop was not able to enforce equilibrium");
00522 if (dt<dtmin)
00523 {
00524 if (Mespr==1) fprintf (stderr,"\n\n time increment is less than minimum time increment");
00525 if (Mespr==1) fprintf (stderr,"\n computation failed (file %s, line %d)\n",__FILE__,__LINE__);
00526
00527
00528
00529 break;
00530 }
00531 }
00532 else
00533 {
00534
00535 if ((Mespr==1) && (j >= 0)) fprintf (stdout,"\n\n Last inner iteration step number j = %ld, norf = %e \n\n",j,norf);
00536
00537
00538
00539 copyv(r, lhsb, n);
00540
00541
00542
00543 copyv(f, fp, n);
00544
00545
00546
00547
00548
00549
00550
00551 Mt->save_nodval(lcid);
00552
00553 Mt->save_nodforce(f);
00554
00555
00556 Gtm->update_auxinf();
00557 prev_timem = Mp->time;
00558
00559 dtdef = Mp->timecon.actualforwtimeincr ();
00560 if (nsts==2)
00561 {
00562 dt*=2.0;
00563 nsts=0;
00564 if (Mespr==1) fprintf (stdout,"\n\n time increment is enlarged because no inner loop was neccessary in previous 3 steps");
00565 }
00566 if (dt>dtdef)
00567 dt=dtdef;
00568
00569
00570
00571
00572
00573
00574
00575
00576 Mm->updateipval();
00577 compute_req_val (lcid);
00578 print_step(lcid, i, Mp->time, f);
00579 print_flush();
00580
00581 Mb->comp_sum_react();
00582 Mb->comp_sum_pdreact();
00583 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], rx=%le, ry=%le, rz=%le",
00584 i, Mp->time, Mp->time/86400, Mb->reactsumcomp[0], Mb->reactsumcomp[1], Mb->reactsumcomp[2]);
00585 fprintf(Out, "\n# i=%ld, time=%le[s], time=%le[d], pdx=%le, pdy=%le, pdz=%le",
00586 i, Mp->time, Mp->time/86400, Mb->pd_reactsumcomp[0], Mb->pd_reactsumcomp[1], Mb->pd_reactsumcomp[2]);
00587 fflush(Out);
00588
00589 if ((Mp->timecon.isitimptime()==1) && (Mp->hdbcont.hdbtype))
00590 {
00591 if (Mespr==1)
00592 fprintf (stdout,"\n Creating backup file\n");
00593 solver_save (r,fp,i,Mp->time,dt,&Mp->timecon,n);
00594 }
00595 }
00596 }while(Mp->time<end_time);
00597
00598 print_close ();
00599
00600 delete [] fi;
00601 delete [] fb;
00602 delete [] fp;
00603 delete [] dr;
00604 delete [] fl;
00605 delete [] lhsb;
00606 delete [] ifn;
00607 }