00001 #include <math.h>
00002 #include <string.h>
00003 #include "stacktrace.h"
00004 #include "mechmat.h"
00005 #include "global.h"
00006 #include "globmat.h"
00007 #include "elemhead.h"
00008 #include "vecttens.h"
00009 #include "intpoints.h"
00010 #include "element.h"
00011 #include "node.h"
00012 #include "elastisomat.h"
00013 #include "elastortomat.h"
00014 #include "elastgmat3d.h"
00015 #include "splas1d.h"
00016 #include "j2flow.h"
00017 #include "j2flow2.h"
00018 #include "microM4.h"
00019 #include "microSIM.h"
00020 #include "microfiber.h"
00021 #include "mohrc.h"
00022 #include "mohrcparab.h"
00023 #include "boermat.h"
00024 #include "drprag.h"
00025 #include "doubdp.h"
00026 #include "drprag2.h"
00027 #include "camclay.h"
00028 #include "camclaycoup.h"
00029 #include "shefplast.h"
00030 #include "hissplas.h"
00031 #include "glasgmech.h"
00032 #include "glasgowdam.h"
00033 #include "creepdam.h"
00034 #include "timeswmat.h"
00035 #include "effstress.h"
00036 #include "simviscous.h"
00037 #include "lemaitre.h"
00038 #include "ortodam.h"
00039 #include "ortodamrot.h"
00040 #include "anisodam.h"
00041 #include "anisodamrot.h"
00042 #include "scaldam.h"
00043 #include "fixortodam.h"
00044 #include "scaldamcc.h"
00045 #include "aci78.h"
00046 #include "cebfip78.h"
00047 #include "graphmat.h"
00048 #include "varelastisomat.h"
00049 #include "geoelast.h"
00050 #include "creep.h"
00051 #include "creep_b3.h"
00052 #include "creep_rspec.h"
00053 #include "creepb.h"
00054 #include "creep_dpl.h"
00055 #include "creep_effym.h"
00056 #include "shrinkmat.h"
00057 #include "consol.h"
00058 #include "winpast.h"
00059 #include "nonlocmicroM4.h"
00060 #include "therisomat.h"
00061 #include "nonlocplast.h"
00062 #include "nonlocdamg.h"
00063 #include "damplast.h"
00064 #include "tensor.h"
00065 #include "visplast.h"
00066 #include "chen.h"
00067 #include "lenjonesmat.h"
00068 #include "contactmat.h"
00069 #include "cebfipcontactmat.h"
00070 #include "elemswitch.h"
00071 #include "relaxeuroc.h"
00072 #include "therisomat.h"
00073 #include "shrinkmat.h"
00074 #include "hypoplast.h"
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 mechmat::mechmat ()
00085 {
00086
00087 nmt=0;
00088
00089 tnip=0;
00090
00091
00092 plast=0;
00093
00094
00095 max_nncompo=0;
00096
00097 max_encompo=0;
00098
00099 max_ncompstrn=0;
00100
00101 max_ncompstre=0;
00102
00103
00104 ip=NULL;
00105
00106 ic=NULL;
00107
00108 elip = NULL;
00109
00110
00111 eliso=NULL; elgm3d=NULL;
00112 spl1d=NULL; j2f=NULL; j2f2=NULL; mcoul=NULL; mcpar=NULL;
00113 boerm=NULL; drprm=NULL; ddpm=NULL; drprm2=NULL; cclay=NULL; cclayc=NULL; shpl=NULL; hisspl=NULL;
00114 glasgmat=NULL; glasgdam=NULL; crdam=NULL; tswmat=NULL; effstr=NULL;
00115 mpM4=NULL; mpSIM=NULL; mpfib=NULL; nmpM4=NULL;
00116 svis=NULL; svipl=NULL; lmtr=NULL;
00117 scdam=NULL; anidam=NULL; anidamrot=NULL; ortdam=NULL;
00118 ortdam=NULL; ortdamrot=NULL;
00119 aci78mod=NULL; cebfip78mod=NULL;
00120 grmat=NULL; geoel=NULL; veliso=NULL;
00121 crb3=NULL;
00122 crrs=NULL;
00123 crdpl=NULL;
00124 crbaz=NULL; csol=NULL; wpast=NULL;
00125 tidilat=NULL; relaxec=NULL;
00126 nlplast=NULL; nldamg=NULL;
00127
00128 dampl = NULL; visplas = NULL;
00129 chplast = NULL;
00130 lenjon = NULL;
00131 conmat=NULL;
00132 cebfipconmat=NULL;
00133
00134
00135 shmat=NULL;
00136
00137
00138 ipv = NULL;
00139
00140
00141 eigstrid = NULL;
00142
00143 eigstrains = NULL;
00144
00145 eigstresses = NULL;
00146
00147 tempstrains = NULL;
00148
00149 nonmechq = NULL;
00150 nnmq = 0;
00151 nmqo = NULL;
00152 for (long i=0; i<tnknmq; i++)
00153 nmqid[i] = -1;
00154
00155 mtype = NULL;
00156 numtype = NULL;
00157 }
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 mechmat::~mechmat ()
00168 {
00169 long i;
00170
00171 delete [] ip;
00172 delete [] elip;
00173 delete [] eliso; delete [] elgm3d;
00174
00175 delete [] spl1d; delete [] j2f; delete [] j2f2;
00176 delete [] mcoul; delete [] mcpar;
00177
00178 delete [] boerm; delete [] drprm; delete [] ddpm; delete [] drprm2; delete [] cclay; delete [] shpl; delete [] hisspl;
00179
00180 delete [] glasgmat; delete [] glasgdam; delete [] crdam; delete [] tswmat; delete [] effstr;
00181
00182 delete [] mpM4; delete [] mpSIM; delete [] mpfib; delete [] nmpM4;
00183
00184 delete [] svis;
00185 delete [] lmtr;
00186
00187 delete [] scdam; delete [] anidam; delete [] anidamrot;
00188 delete [] ortdam; delete [] ortdamrot;
00189
00190 delete [] aci78mod; delete [] cebfip78mod;
00191 delete [] grmat; delete [] geoel; delete [] veliso;
00192 delete [] crb3; delete [] crrs; delete [] crdpl;
00193 delete [] crbaz; delete [] csol; delete [] wpast;
00194 delete [] tidilat; delete [] relaxec;
00195 delete [] nlplast; delete [] nldamg;
00196 delete [] dampl; delete [] visplas;
00197 delete [] chplast;
00198 delete [] lenjon; delete [] conmat; delete [] cebfipconmat;
00199
00200
00201 delete [] shmat;
00202
00203 if (ic != NULL)
00204 {
00205 for (i = 0; i < tnip; i++)
00206 delete [] ic[i];
00207 }
00208 delete [] ic;
00209
00210 delete [] ipv;
00211
00212 if (nonmechq != NULL){
00213 for (i=0;i<nnmq;i++){
00214 delete [] nonmechq[i];
00215 }
00216 delete [] nonmechq;
00217 }
00218
00219 delete [] nmqo;
00220
00221 if (eigstrains != NULL){
00222 for (i=0;i<tnip;i++){
00223 delete [] eigstrains[i];
00224 }
00225 delete [] eigstrains;
00226 }
00227
00228 if (eigstrid != NULL){
00229 for (i=0;i<tnip;i++){
00230 delete [] eigstrid[i];
00231 }
00232 delete [] eigstrid;
00233 }
00234
00235 if (eigstresses != NULL){
00236 for (i=0;i<tnip;i++){
00237 delete [] eigstresses[i];
00238 }
00239 delete [] eigstresses;
00240 }
00241 if (tempstrains != NULL){
00242 for (i=0;i<tnip;i++){
00243 delete [] tempstrains[i];
00244 }
00245 delete [] tempstrains;
00246 }
00247
00248 delete [] mtype;
00249 delete [] numtype;
00250 }
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 long mechmat::intpnum (void)
00263 {
00264 long i,j,k,n,nb;
00265
00266 n=0;
00267 for (i=0;i<Mt->ne;i++){
00268 nb=Mt->give_nb (i);
00269 Mt->elements[i].ipp = new long* [nb];
00270 for (j=0;j<nb;j++){
00271 Mt->elements[i].ipp[j] = new long [nb];
00272 }
00273 for (j=0;j<nb;j++){
00274 for (k=0;k<nb;k++){
00275 Mt->elements[i].ipp[j][k]=n;
00276 n+=Mt->give_nip (i,j,k);
00277 }
00278 }
00279 }
00280 return n;
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 void mechmat::intpointalloc ()
00293 {
00294 ip = new intpoints [tnip];
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 void mechmat::alloc_ic ()
00307 {
00308 long i,n;
00309
00310 ic = new double* [tnip];
00311 for (i=0;i<tnip;i++){
00312 n = ip[i].ncompstr;
00313 ic[i] = new double[3+n];
00314 }
00315 }
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 void mechmat::alloc_nonmechq(long n)
00331 {
00332 long i, j;
00333
00334 nnmq = n;
00335 nonmechq = new double*[nnmq];
00336
00337 for (i=0; i<nnmq; i++)
00338 {
00339 nonmechq[i] = new double[tnip];
00340
00341 memset(nonmechq[i], 0, sizeof(*nonmechq[i])*tnip);
00342 }
00343 if (nmqo == NULL)
00344 {
00345 nmqo = new nonmechquant[nnmq];
00346 j=0;
00347 for (i=0; i<tnknmq; i++)
00348 {
00349 if (nmqid[i] > -1)
00350 {
00351 nmqo[j] = nonmechquant(i+1);
00352 j++;
00353 }
00354 }
00355 if (j != nnmq)
00356 {
00357 print_err("indices of used non-mechanical qunatities must be initialized before allocation", __FILE__, __LINE__, __func__);
00358 abort();
00359 }
00360 }
00361 return;
00362 }
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 void mechmat::readip (XFILE *in)
00379 {
00380 long i;
00381 for (i=0;i<tnip;i++){
00382 ip[i].read (in);
00383 }
00384 }
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396 void mechmat::intpointinit ()
00397 {
00398 long i, j, ii, jj, k, nb, nip, ipp;
00399
00400
00401 elip = new long[tnip];
00402
00403 for (i=0;i<Mt->ne;i++){
00404 nb=Mt->give_nb (i);
00405 for (ii=0;ii<nb;ii++){
00406 for (jj=0;jj<nb;jj++){
00407 nip=Mt->give_nip(i,ii,jj);
00408 ipp=Mt->elements[i].ipp[ii][jj];
00409 for (j=0;j<nip;j++){
00410 if (elip != NULL)
00411 elip[ipp] = i;
00412 ip[ipp].nm = Mt->elements[i].nm;
00413 ip[ipp].tm = new mattype[ip[ipp].nm];
00414 ip[ipp].idm = new long[ip[ipp].nm];
00415 for (k = 0; k < ip[ipp].nm; k++)
00416 {
00417 ip[ipp].tm[k] = Mt->elements[i].tm[k];
00418 ip[ipp].idm[k] = Mt->elements[i].idm[k];
00419 switch (ip[ipp].tm[k])
00420 {
00421 case therisodilat:
00422 ip[ipp].hmt ^= 1;
00423 break;
00424 case nonlocplastmat:
00425 ip[ipp].hmt ^= 2;
00426 break;
00427 case nonlocdamgmat:
00428 ip[ipp].hmt ^= 2;
00429 break;
00430 case nonlocalmicroM4:
00431 ip[ipp].hmt ^= 2;
00432 break;
00433 case relaxationeuro:
00434 ip[ipp].hmt ^= 4;
00435 break;
00436 case effective_stress:
00437 ip[ipp].hmt ^= 8;
00438 break;
00439 default:
00440 break;
00441 }
00442 }
00443 ipp++;
00444 }
00445 }
00446 }
00447 }
00448 }
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460 void mechmat::init_ip_1 (void)
00461 {
00462 long i,j,ne,ii,jj,ipp,nb,nip,ncomp;
00463 strastrestate ssst;
00464
00465 ne=Mt->ne;
00466
00467 for (i=0;i<ne;i++){
00468 nb=Mt->give_nb (i);
00469 for (ii=0;ii<nb;ii++){
00470
00471 ncomp=Mt->give_ncomp (i,ii);
00472
00473 ssst=Mt->give_ssst (i,ii);
00474 for (jj=0;jj<nb;jj++){
00475
00476 ipp=Mt->elements[i].ipp[ii][jj];
00477
00478 nip=Mt->give_nip (i,ii,jj);
00479 for (j=0;j<nip;j++){
00480 ip[ipp].ncompstr=ncomp;
00481 ip[ipp].ssst=ssst;
00482
00483 if ((ssst == planestrain) || (ssst == planestress))
00484 ip[ipp].ncompstr=4;
00485
00486 ipp++;
00487 }
00488 }
00489 }
00490 }
00491 }
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507 void mechmat::init_ip_2 (void)
00508 {
00509 long i,n,ii,nlmid,ncompnl;
00510
00511
00512
00513
00514 if ((Mp->temperature == 1) || (Mp->pore_press == 1))
00515 {
00516
00517 n = search_reqnmq(nmqo);
00518 if (n)
00519 {
00520 for (i=0; i<n; i++)
00521 nmqid[nmqo[i]-1] = i;
00522
00523 alloc_nonmechq(n);
00524 }
00525 }
00526 else
00527 {
00528
00529 }
00530
00531
00532 for (i=0;i<tnip;i++)
00533 {
00534 ncompnl = 0;
00535 ii = ip[i].hmt & 2;
00536 if (ii > 0)
00537
00538 {
00539 nlmid = givenonlocid(i);
00540 ncompnl = give_num_averq(i, nlmid);
00541 }
00542
00543 ip[i].alloc(Mb->nlc,i,ncompnl);
00544 }
00545 }
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559 void mechmat::read (XFILE *in)
00560 {
00561
00562
00563 readmatchar (in);
00564
00565
00566 tnip = intpnum ();
00567 if (Mespr==1){
00568 fprintf (stdout,"\n number of integration points %ld",tnip);
00569 }
00570
00571
00572 intpointalloc ();
00573
00574
00575
00576 intpointinit ();
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 }
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602 void mechmat::readmatchar (XFILE *in)
00603 {
00604 long i;
00605
00606 xfscanf (in, "%k%ld", "num_mat_types", &nmt);
00607
00608 numtype = new long [nmt];
00609 mtype = new mattype [nmt];
00610
00611 if (nmt<1)
00612 print_err("wrong number of material types",__FILE__,__LINE__,__func__);
00613
00614 if (Mespr==1) fprintf (stdout,"\n number of different types of materials %ld",nmt);
00615
00616 for (i=0;i<nmt;i++)
00617 {
00618 xfscanf (in,"%k%m %k%ld", "mattype", &mattype_kwdset, &mtype[i], "num_inst", &numtype[i]);
00619 if (numtype[i]<1)
00620 print_err("wrong number of material characteristics",__FILE__,__LINE__,__func__);
00621
00622 readmattype(in, mtype[i], numtype[i]);
00623 }
00624 }
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639 void mechmat::readmattype(XFILE *in, mattype mtype, long numt)
00640 {
00641 long j,k;
00642
00643 switch (mtype)
00644 {
00645 case elisomat:{
00646 if (Mespr==1) fprintf (stdout,"\n number of elastic isotropic materials %ld",numt);
00647 eliso = new elastisomat [numt];
00648 for (j=0;j<numt;j++){
00649 k=numt+1;
00650 xfscanf (in,"%ld",&k);
00651 if (k>numt || k<1)
00652 print_err("wrong number of elastic isotropic material",__FILE__,__LINE__,__func__);
00653 eliso[k-1].read (in);
00654 }
00655 break;
00656 }
00657 case elortomat:{
00658 if (Mespr==1) fprintf (stdout,"\n number of elastic orthotropic materials %ld",numt);
00659 elorto = new elastortomat [numt];
00660 for (j=0;j<numt;j++){
00661 k=numt+1;
00662 xfscanf (in,"%ld",&k);
00663 if (k>numt || k<1)
00664 print_err("wrong number of elastic orthtropic material",__FILE__,__LINE__,__func__);
00665 elorto[k-1].read (in);
00666 }
00667 break;
00668 }
00669 case elgmat3d:{
00670 if (Mespr==1) fprintf (stdout,"\n number of elastic anisotropic materials %ld",numt);
00671 elgm3d = new elastgmat3d [numt];
00672 for (j=0;j<numt;j++){
00673 k=numt+1;
00674 xfscanf (in,"%ld",&k);
00675 if (k>numt || k<1)
00676 print_err("wrong number of elastic anisotropic material",__FILE__,__LINE__,__func__);
00677 elgm3d[k-1].read (in);
00678 }
00679 break;
00680 }
00681
00682 case simplas1d:{
00683 if (Mespr==1) fprintf (stdout,"\n number of simple plastic 1D materials %ld",numt);
00684 spl1d = new splas1d [numt];
00685 for (j=0;j<numt;j++){
00686 k=numt+1;
00687 xfscanf (in,"%ld",&k);
00688 if (k>numt || k<1)
00689 print_err("wrong number of simple plastic 1D material",__FILE__,__LINE__,__func__);
00690 spl1d[k-1].read (in);
00691 }
00692 break;
00693 }
00694 case jflow:{
00695 if (Mespr==1) fprintf (stdout,"\n number of J2 flow materials %ld",numt);
00696 j2f = new j2flow [numt];
00697 for (j=0;j<numt;j++){
00698 k=numt+1;
00699 xfscanf (in,"%ld",&k);
00700 if (k>numt || k<1)
00701 print_err("wrong number of J2 flow material",__FILE__,__LINE__,__func__);
00702 j2f[k-1].read (in);
00703 }
00704 break;
00705 }
00706 case jflow2:{
00707 if (Mespr==1) fprintf (stdout,"\n number of J2 flow materials %ld",numt);
00708 j2f2 = new j2flow2 [numt];
00709 for (j=0;j<numt;j++){
00710 k=numt+1;
00711 xfscanf (in,"%ld",&k);
00712 if (k>numt || k<1)
00713 print_err("wrong number of J2 flow material",__FILE__,__LINE__,__func__);
00714 j2f2[k-1].read (in);
00715 }
00716 break;
00717 }
00718 case microplaneM4:{
00719 if (Mespr==1) fprintf (stdout,"\n number of microplane materials %ld",numt);
00720 mpM4 = new microM4 [numt];
00721 for (j=0;j<numt;j++){
00722 k=numt+1;
00723 xfscanf (in,"%ld",&k);
00724 if (k>numt || k<1)
00725 print_err("wrong number of microplane material",__FILE__,__LINE__,__func__);
00726 mpM4[k-1].read (in);
00727 }
00728 break;
00729 }
00730 case nonlocalmicroM4:{
00731 if (Mespr==1) fprintf (stdout,"\n number of nonlocal microplane materials %ld",numt);
00732 nmpM4 = new nonlocmicroM4 [numt];
00733 for (j=0;j<numt;j++){
00734 k=numt+1;
00735 xfscanf (in,"%ld",&k);
00736 if (k>numt || k<1)
00737 print_err("wrong number of nonlocal microplane material",__FILE__,__LINE__,__func__);
00738 nmpM4[k-1].read (in);
00739 }
00740 break;
00741 }
00742 case microsimp:{
00743 if (Mespr==1) fprintf (stdout,"\n number of microplane materials %ld",numt);
00744 mpSIM = new microSIM [numt];
00745 for (j=0;j<numt;j++){
00746 k=numt+1;
00747 xfscanf (in,"%ld",&k);
00748 if (k>numt || k<1)
00749 print_err("wrong number of microplane material",__FILE__,__LINE__,__func__);
00750 mpSIM[k-1].read (in);
00751 }
00752 break;
00753 }
00754 case microfibro:{
00755 if (Mespr==1) fprintf (stdout,"\n number of microplane materials %ld",numt);
00756 mpfib = new microfiber [numt];
00757 for (j=0;j<numt;j++){
00758 k=numt+1;
00759 xfscanf (in,"%ld",&k);
00760 if (k>numt || k<1)
00761 print_err("wrong number of microplane material",__FILE__,__LINE__,__func__);
00762 mpfib[k-1].read (in);
00763 }
00764 break;
00765 }
00766 case mohrcoul:{
00767 if (Mespr==1) fprintf (stdout,"\n number of Mohr-Coulomb materials %ld",numt);
00768 mcoul = new mohrcoulomb [numt];
00769 for (j=0;j<numt;j++){
00770 k=numt+1;
00771 xfscanf (in,"%ld",&k);
00772 if (k>numt || k<1)
00773 print_err("wrong number of Mohr-Coulomb material",__FILE__,__LINE__,__func__);
00774 mcoul[k-1].read (in);
00775 }
00776 break;
00777 }
00778 case mohrcoulparab:{
00779 if (Mespr==1) fprintf (stdout,"\n number of parabolic Mohr-Coulomb materials %ld",numt);
00780 mcpar = new mohrcoulombpar [numt];
00781 for (j=0;j<numt;j++){
00782 k=numt+1;
00783 xfscanf (in,"%ld",&k);
00784 if (k>numt || k<1)
00785 print_err("wrong number of parabolic Mohr-Coulomb material",__FILE__,__LINE__,__func__);
00786 mcpar[k-1].read (in);
00787 }
00788 break;
00789 }
00790 case boermaterial:{
00791 if (Mespr==1) fprintf (stdout,"\n number of Boer materials %ld",numt);
00792 boerm = new boermat [numt];
00793 for (j=0;j<numt;j++){
00794 k=numt+1;
00795 xfscanf (in,"%ld",&k);
00796 if (k>numt || k<1)
00797 print_err("wrong number of Boer material",__FILE__,__LINE__,__func__);
00798 boerm[k-1].read (in);
00799 }
00800 break;
00801 }
00802 case druckerprager:{
00803 if (Mespr==1) fprintf (stdout,"\n number of Drucker-Prager materials %ld",numt);
00804 drprm = new drprag [numt];
00805 for (j=0;j<numt;j++){
00806 k=numt+1;
00807 xfscanf (in,"%ld",&k);
00808 if (k>numt || k<1)
00809 print_err("wrong number of Drucker-Prager material",__FILE__,__LINE__,__func__);
00810 drprm[k-1].read (in);
00811 }
00812 break;
00813 }
00814 case doubledrprager:{
00815 if (Mespr==1) fprintf (stdout,"\n number of Drucker-Prager materials %ld",numt);
00816 ddpm = new doubdp [numt];
00817 for (j=0;j<numt;j++){
00818 k=numt+1;
00819 xfscanf (in,"%ld",&k);
00820 if (k>numt || k<1)
00821 print_err("wrong number of Drucker-Prager material",__FILE__,__LINE__,__func__);
00822 ddpm[k-1].read (in);
00823 }
00824 break;
00825 }
00826 case druckerprager2:{
00827 if (Mespr==1) fprintf (stdout,"\n number of Drucker-Prager2 materials %ld",numt);
00828 drprm2 = new drprag2 [numt];
00829 for (j=0;j<numt;j++){
00830 k=numt+1;
00831 xfscanf (in,"%ld",&k);
00832 if (k>numt || k<1)
00833 print_err("wrong number of Drucker-Prager material",__FILE__,__LINE__,__func__);
00834 drprm2[k-1].read (in);
00835 }
00836 break;
00837 }
00838 case modcamclaymat:{
00839 if (Mespr==1) fprintf (stdout,"\n number of modified cam-clay materials %ld",numt);
00840 cclay = new camclay [numt];
00841 for (j=0;j<numt;j++){
00842 k=numt+1;
00843 xfscanf (in,"%ld",&k);
00844 if (k>numt || k<1)
00845 print_err("wrong number of modified cam-clay material",__FILE__,__LINE__,__func__);
00846 cclay[k-1].read (in);
00847 }
00848 break;
00849 }
00850 case modcamclaycoupmat:{
00851 if (Mespr==1) fprintf (stdout,"\n number of modified cam-clay materials for coupling %ld",numt);
00852 cclayc = new camclaycoup [numt];
00853 for (j=0;j<numt;j++){
00854 k=numt+1;
00855 xfscanf (in,"%ld",&k);
00856 if (k>numt || k<1)
00857 print_err("wrong number of modified cam-clay material for coupling",__FILE__,__LINE__,__func__);
00858 cclayc[k-1].read (in);
00859 }
00860 break;
00861 }
00862 case hypoplastmat:{
00863 if (Mespr==1) fprintf (stdout,"\n number of modified hypoplast materials for coupling %ld",numt);
00864 hypopl = new hypoplast [numt];
00865 for (j=0;j<numt;j++){
00866 k=numt+1;
00867 xfscanf (in,"%ld",&k);
00868 if (k>numt || k<1)
00869 print_err("wrong number of hypoplast material for coupling",__FILE__,__LINE__,__func__);
00870 hypopl[k-1].read (in);
00871 }
00872 break;
00873 }
00874 case shefpl:{
00875 if (Mespr==1) fprintf (stdout,"\n number of Sheffield plasticity materials %ld",numt);
00876 shpl = new shefplast [numt];
00877 for (j=0;j<numt;j++){
00878 k=numt+1;
00879 xfscanf (in,"%ld",&k);
00880 if (k>numt || k<1)
00881 print_err("wrong number of Sheffiled material",__FILE__,__LINE__,__func__);
00882 shpl[k-1].read (in);
00883 }
00884 break;
00885 }
00886 case hissplasmat:{
00887 if (Mespr==1) fprintf (stdout,"\n number of HISS plasticity materials %ld",numt);
00888 hisspl = new hissplas [numt];
00889 for (j=0;j<numt;j++){
00890 k=numt+1;
00891 xfscanf (in,"%ld",&k);
00892 if (k>numt || k<1)
00893 print_err("wrong number of HISS plasticity material",__FILE__,__LINE__,__func__);
00894 hisspl[k-1].read (in);
00895 }
00896 break;
00897 }
00898 case chenplast:{
00899 if (Mespr==1) fprintf (stdout,"\n number of Chen plasticity materials %ld",numt);
00900 chplast = new chen [numt];
00901 for (j=0;j<numt;j++){
00902 k=numt+1;
00903 xfscanf (in,"%ld",&k);
00904 if (k>numt || k<1)
00905 print_err("wrong number of Chen material",__FILE__,__LINE__,__func__);
00906 chplast[k-1].read (in);
00907 }
00908 break;
00909 }
00910 case glasgowmechmat:{
00911 if (Mespr==1) fprintf (stdout,"\n number of Glasgow plasticity materials %ld",numt);
00912 glasgmat = new glasgmech [numt];
00913 for (j=0;j<numt;j++){
00914 k=numt+1;
00915 xfscanf (in,"%ld",&k);
00916 if (k>numt || k<1)
00917 print_err("wrong number of Glasgow material",__FILE__,__LINE__,__func__);
00918 glasgmat[k-1].read (in);
00919 }
00920 break;
00921 }
00922 case glasgowdamage:{
00923 if (Mespr==1) fprintf (stdout,"\n number of Glasgow damage materials %ld",numt);
00924 glasgdam = new glasgowdam [numt];
00925 for (j=0;j<numt;j++){
00926 k=numt+1;
00927 xfscanf (in,"%ld",&k);
00928 if (k>numt || k<1)
00929 print_err("wrong number of Glasgow damage material",__FILE__,__LINE__,__func__);
00930 glasgdam[k-1].read (in);
00931 }
00932 break;
00933 }
00934 case creep_damage:{
00935 if (Mespr==1) fprintf (stdout,"\n number of artificial material creep-damage %ld",numt);
00936 crdam = new creepdam [numt];
00937 for (j=0;j<numt;j++){
00938 k=numt+1;
00939 xfscanf (in,"%ld",&k);
00940 if (k>numt || k<1)
00941 print_err("wrong number of artificial material creep-damge",__FILE__,__LINE__,__func__);
00942 }
00943 break;
00944 }
00945 case time_switchmat:{
00946 if (Mespr==1) fprintf (stdout,"\n number of artificial time-switched material %ld",numt);
00947 tswmat = new timeswmat [numt];
00948 for (j=0;j<numt;j++){
00949 k=numt+1;
00950 xfscanf (in,"%ld",&k);
00951 if (k>numt || k<1)
00952 print_err("wrong number of artificial time-switched material",__FILE__,__LINE__,__func__);
00953 }
00954 break;
00955 }
00956 case effective_stress:{
00957 if (Mespr==1) fprintf (stdout,"\n number of artificial material for effective stresses %ld",numt);
00958 effstr = new effstress [numt];
00959 for (j=0;j<numt;j++){
00960 k=numt+1;
00961 xfscanf (in,"%ld",&k);
00962 if (k>numt || k<1)
00963 print_err("wrong number of artificial material for effective stresses",__FILE__,__LINE__,__func__);
00964 effstr[k-1].read (in);
00965 }
00966 break;
00967 }
00968 case simvisc:{
00969 if (Mespr==1) fprintf (stdout,"\n number of simple viscous materials %ld",numt);
00970 svis = new simviscous [numt];
00971 for (j=0;j<numt;j++){
00972 k=numt+1;
00973 xfscanf (in,"%ld",&k);
00974 if (k>numt || k<1)
00975 print_err("wrong number of simple viscous material in function",__FILE__,__LINE__,__func__);
00976 svis[k-1].read (in);
00977 }
00978 break;
00979 }
00980 case lemaitr:{
00981 if (Mespr==1) fprintf (stdout,"\n number of Lemaitre visco materials %ld",numt);
00982 lmtr = new lemaitre [numt];
00983 for (j=0;j<numt;j++){
00984 k=numt+1;
00985 xfscanf (in,"%ld",&k);
00986 if (k>numt || k<1)
00987 print_err("wrong number of Lemaitre visco material",__FILE__,__LINE__,__func__);
00988 lmtr[k-1].read (in);
00989 }
00990 break;
00991 }
00992 case scaldamage:{
00993 if (Mespr==1) fprintf (stdout,"\n number of simple damage materials %ld",numt);
00994 scdam = new scaldam [numt];
00995 for (j=0;j<numt;j++){
00996 k=numt+1;
00997 xfscanf (in,"%ld",&k);
00998 if (k>numt || k<1)
00999 print_err("wrong number of simple damage material",__FILE__,__LINE__,__func__);
01000 scdam[k-1].read (in);
01001 }
01002 break;
01003 }
01004 case fixortodamage:{
01005 if (Mespr==1) fprintf (stdout,"\n number of orthotropic damage materials %ld",numt);
01006 fixodam = new fixortodam [numt];
01007 for (j=0;j<numt;j++){
01008 k=numt+1;
01009 xfscanf (in,"%ld",&k);
01010 if (k>numt || k<1)
01011 print_err("wrong number of orthotropic damage material",__FILE__,__LINE__,__func__);
01012 fixodam[k-1].read (in);
01013 }
01014 break;
01015 }
01016 case scaldamagecc:{
01017 if (Mespr==1) fprintf (stdout,"\n number of simple damage materials with crack closure %ld",numt);
01018 scdamcc = new scaldamcc [numt];
01019 for (j=0;j<numt;j++){
01020 k=numt+1;
01021 xfscanf (in,"%ld",&k);
01022 if (k>numt || k<1)
01023 print_err("wrong number of simple damage material with crack closure",__FILE__,__LINE__,__func__);
01024 scdamcc[k-1].read (in);
01025 }
01026 break;
01027 }
01028 case ortodamage:{
01029 if (Mespr==1) fprintf (stdout,"\n number of ortotropic damage materials %ld",numt);
01030 ortdam = new ortodam [numt];
01031 for (j=0;j<numt;j++){
01032 k=numt+1;
01033 xfscanf (in,"%ld",&k);
01034 if (k>numt || k<1)
01035 print_err("wrong number of ortotropic damage material",__FILE__,__LINE__,__func__);
01036 ortdam[k-1].read (in);
01037 }
01038 break;
01039 }
01040 case ortodamagerot:{
01041 if (Mespr==1) fprintf (stdout,"\n number of orthotropic damage materials with crack rotation %ld",numt);
01042 ortdamrot = new ortodamrot [numt];
01043 for (j=0;j<numt;j++){
01044 k=numt+1;
01045 xfscanf (in,"%ld",&k);
01046 if (k>numt || k<1)
01047 print_err("wrong number of orthotropic damage material with crack rotation",__FILE__,__LINE__,__func__);
01048 ortdamrot[k-1].read (in);
01049 }
01050 break;
01051 }
01052 case anisodamage:{
01053 if (Mespr==1) fprintf (stdout,"\n number of anisotropic damage materials %ld",numt);
01054 anidam = new anisodam [numt];
01055 for (j=0;j<numt;j++){
01056 k=numt+1;
01057 xfscanf (in,"%ld",&k);
01058 if (k>numt || k<1)
01059 print_err("wrong number of anisotropic damage material",__FILE__,__LINE__,__func__);
01060 anidam[k-1].read (in);
01061 }
01062 break;
01063 }
01064 case anisodamagerot:{
01065 if (Mespr==1) fprintf (stdout,"\n number of rotational anisotropic damage materials with crack rotation %ld",numt);
01066 anidamrot = new anisodamrot [numt];
01067 for (j=0;j<numt;j++){
01068 k=numt+1;
01069 xfscanf (in,"%ld",&k);
01070 if (k>numt || k<1)
01071 print_err("wrong number of anisotropic damage material with crack rotation",__FILE__,__LINE__,__func__);
01072 anidamrot[k-1].read (in);
01073 }
01074 break;
01075 }
01076
01077 case aci:{
01078 if (Mespr==1) fprintf (stdout,"\n number of aci78 materials %ld",numt);
01079 aci78mod = new aci78 [numt];
01080 for (j=0;j<numt;j++){
01081 k=numt+1;
01082 xfscanf (in,"%ld",&k);
01083 if (k>numt || k<1)
01084 print_err("wrong number of aci78 material",__FILE__,__LINE__,__func__);
01085 aci78mod[k-1].read (in);
01086 }
01087 break;
01088 }
01089 case cebfip:{
01090 if (Mespr==1) fprintf (stdout,"\n number of cebfip78 materials %ld",numt);
01091 cebfip78mod = new cebfip78 [numt];
01092 for (j=0;j<numt;j++){
01093 k=numt+1;
01094 xfscanf (in,"%ld",&k);
01095 if (k>numt || k<1)
01096 print_err("wrong number of cebfip78 material in function",__FILE__,__LINE__,__func__);
01097 cebfip78mod[k-1].read (in);
01098 }
01099 break;
01100 }
01101 case graphm:{
01102 if (Mespr==1) fprintf (stdout,"\n number of graphic material materials %ld",numt);
01103 grmat = new graphmat [numt];
01104 for (j=0;j<numt;j++){
01105 k=numt+1;
01106 xfscanf (in,"%ld",&k);
01107 if (k>numt || k<1)
01108 print_err("wrong number of graphic material",__FILE__,__LINE__,__func__);
01109 grmat[k-1].read (in);
01110 }
01111 break;
01112 }
01113 case varelisomat:{
01114 if (Mespr==1) fprintf (stdout,"\n number of variable elastic isotropic materials %ld",numt);
01115 veliso = new varelastisomat [numt];
01116 for (j=0;j<numt;j++){
01117 k=numt+1;
01118 xfscanf (in,"%ld",&k);
01119 if (k>numt || k<1)
01120 print_err("wrong number of variable elastic isotropic material",__FILE__,__LINE__,__func__);
01121 veliso[k-1].read (in);
01122 }
01123 break;
01124 }
01125 case geoelast:{
01126 if (Mespr==1) fprintf (stdout,"\n number of geoelast material materials %ld",numt);
01127 geoel = new geoelastmat [numt];
01128 for (j=0;j<numt;j++){
01129 k=numt+1;
01130 xfscanf (in,"%ld",&k);
01131 if (k>numt || k<1)
01132 print_err("wrong number of graphic material",__FILE__,__LINE__,__func__);
01133 geoel[k-1].read (in);
01134 }
01135 break;
01136 }
01137 case creepb3:{
01138 if (Mespr==1) fprintf (stdout,"\n number of b3 creep material materials %ld",numt);
01139 crb3 = new b3mat [numt];
01140 for (j=0;j<numt;j++){
01141 k=numt+1;
01142 xfscanf (in,"%ld",&k);
01143 if (k>numt || k<1)
01144 print_err("wrong number of b3 creep material",__FILE__,__LINE__,__func__);
01145 crb3[k-1].read (in);
01146 }
01147 break;
01148 }
01149 case creeprs:{
01150 if (Mespr==1) fprintf (stdout,"\n number of b3 creep material materials %ld",numt);
01151 crrs = new rspecmat [numt];
01152 for (j=0;j<numt;j++){
01153 k=numt+1;
01154 xfscanf (in,"%ld",&k);
01155 if (k>numt || k<1)
01156 print_err("wrong number of b3 creep material",__FILE__,__LINE__,__func__);
01157 crrs[k-1].read (in);
01158 }
01159 break;
01160 }
01161 case creepdpl:{
01162 if (Mespr==1) fprintf (stdout,"\n number of double-power law creep material materials %ld",numt);
01163 crdpl = new dplmat [numt];
01164 for (j=0;j<numt;j++){
01165 k=numt+1;
01166 xfscanf (in,"%ld",&k);
01167 if (k>numt || k<1)
01168 print_err("wrong number of double-power law creep material",__FILE__,__LINE__,__func__);
01169 crdpl[k-1].read (in);
01170 }
01171 break;
01172 }
01173 case creepbaz:{
01174 if (Mespr==1) fprintf (stdout,"\n number of creep material materials %ld",numt);
01175 crbaz = new creepb [numt];
01176 for (j=0;j<numt;j++){
01177 k=numt+1;
01178 xfscanf (in,"%ld",&k);
01179 if (k>numt || k<1)
01180 print_err("wrong number of creep material in function",__FILE__,__LINE__,__func__);
01181 crbaz[k-1].read (in);
01182 }
01183 break;
01184 }
01185 case creepeffym:{
01186 if (Mespr==1) fprintf (stdout,"\n number of creep_effym materials %ld",numt);
01187 creffym = new creep_effym [numt];
01188 for (j=0;j<numt;j++){
01189 k=numt+1;
01190 xfscanf (in,"%ld",&k);
01191 if (k>numt || k<1)
01192 print_err ("wrong number of creep_effym material in function",__FILE__,__LINE__,__func__);
01193 }
01194 break;
01195 }
01196 case consolidation:{
01197 if (Mespr==1) fprintf (stdout,"\n number of consolidation material materials %ld",numt);
01198 csol = new consol [numt];
01199 for (j=0;j<numt;j++){
01200 k=numt+1;
01201 xfscanf (in,"%ld",&k);
01202 if (k>numt || k<1)
01203 print_err("wrong number of consolidation material model",__FILE__,__LINE__,__func__);
01204 csol[k-1].read (in);
01205 }
01206 break;
01207 }
01208 case winklerpasternak:{
01209 if (Mespr==1) fprintf (stdout,"\n number of Winkler-Pasternak material materials %ld",numt);
01210 wpast = new winpast [numt];
01211 for (j=0;j<numt;j++){
01212 k=numt+1;
01213 xfscanf (in,"%ld",&k);
01214 if (k>numt || k<1)
01215 print_err("wrong number of Winkler-Pasternak material model",__FILE__,__LINE__,__func__);
01216 wpast[k-1].read (in);
01217 }
01218 break;
01219 }
01220 case therisodilat:{
01221 if (Mespr==1) fprintf (stdout,"\n number of isotropic thermal dilatancy materials %ld",numt);
01222 tidilat = new therisomat [numt];
01223 for (j=0;j<numt;j++){
01224 k=numt+1;
01225 xfscanf (in,"%ld",&k);
01226 if (k>numt || k<1)
01227 print_err("wrong number of isotropic thermal dilatancy material",__FILE__,__LINE__,__func__);
01228 tidilat[k-1].read (in);
01229 }
01230 break;
01231 }
01232 case relaxationeuro:{
01233 if (Mespr==1) fprintf (stdout,"\n number of models of stress relaxation based on Eurocode2 %ld",numt);
01234 relaxec = new relaxeuroc [numt];
01235 for (j=0;j<numt;j++){
01236 k=numt+1;
01237 xfscanf (in,"%ld",&k);
01238 if (k>numt || k<1)
01239 print_err("wrong number of models of stress relaxation based on Eurocode2",__FILE__,__LINE__,__func__);
01240 relaxec[k-1].read (in);
01241 }
01242 break;
01243 }
01244
01245 case nonlocplastmat:{
01246 if (Mespr==1) fprintf (stdout,"\n number of nonlocal plasticity materials %ld",numt);
01247 nlplast = new nonlocplast [numt];
01248 for (j=0;j<numt;j++){
01249 k=numt+1;
01250 xfscanf (in,"%ld",&k);
01251 if (k>numt || k<1)
01252 print_err("wrong number of nonlocal plasticity material",__FILE__,__LINE__,__func__);
01253 nlplast[k-1].read (in);
01254 }
01255
01256 Mp->matmodel = nonlocal;
01257 break;
01258 }
01259 case nonlocdamgmat:{
01260 if (Mespr==1) fprintf (stdout,"\n number of nonlocal damage materials %ld",numt);
01261 nldamg = new nonlocdamg [numt];
01262 for (j=0;j<numt;j++){
01263 k=numt+1;
01264 xfscanf (in,"%ld",&k);
01265 if (k>numt || k<1)
01266 print_err("wrong number of nonlocal damage material", __FILE__,__LINE__,__func__);
01267 nldamg[k-1].read (in);
01268 }
01269
01270 Mp->matmodel = nonlocal;
01271 break;
01272 }
01273
01274 case damage_plasticity:{
01275 if (Mespr==1) fprintf (stdout,"\n number of artificial material damage-plasticity %ld",numt);
01276 dampl = new damplast [numt];
01277 for (j=0;j<numt;j++){
01278 k=numt+1;
01279 xfscanf (in,"%ld",&k);
01280 if (k>numt || k<1)
01281 print_err("wrong number of artificial material damge-plasticity",__FILE__,__LINE__,__func__);
01282 }
01283 break;
01284 }
01285
01286 case viscoplasticity:{
01287 if (Mespr==1) fprintf (stdout,"\n number of artificial material of viscoplasticity %ld",numt);
01288 visplas = new visplast [numt];
01289 for (j=0;j<numt;j++){
01290 k=numt+1;
01291 xfscanf (in,"%ld",&k);
01292 if (k>numt || k<1)
01293 print_err("wrong number of artificial material of viscoplasticity",__FILE__,__LINE__,__func__);
01294 }
01295 break;
01296 }
01297
01298 case lenjonespot:{
01299 if (Mespr==1) fprintf (stdout,"\n number of Lennard-Jones interatomic potentials %ld",numt);
01300 lenjon = new lenjonesmat [numt];
01301 for (j=0;j<numt;j++){
01302 k=numt+1;
01303 xfscanf (in,"%ld",&k);
01304 if (k>numt || k<1)
01305 print_err("wrong number of Lennard-Jones interatomic potentials",__FILE__,__LINE__,__func__);
01306 lenjon[k-1].read (in);
01307 }
01308 break;
01309 }
01310
01311 case contmat:{
01312 if (Mespr==1) fprintf (stdout,"\n number of materials for contact problems %ld",numt);
01313 conmat = new contactmat [numt];
01314 for (j=0;j<numt;j++){
01315 k=numt+1;
01316 xfscanf (in,"%ld",&k);
01317 if (k>numt || k<1)
01318 print_err("wrong number of materials for contact problems",__FILE__,__LINE__,__func__);
01319 conmat[k-1].read (in);
01320 }
01321 break;
01322 }
01323
01324 case cebfipcontmat:{
01325 if (Mespr==1) fprintf (stdout,"\n number of materials for contact problems %ld",numt);
01326 cebfipconmat = new cebfipcontactmat [numt];
01327 for (j=0;j<numt;j++){
01328 k=numt+1;
01329 xfscanf (in,"%ld",&k);
01330 if (k>numt || k<1)
01331 print_err("wrong number of materials for contact problems",__FILE__,__LINE__,__func__);
01332 cebfipconmat[k-1].read (in);
01333 }
01334 break;
01335 }
01336
01337 case shrinkagemat:{
01338 if (Mespr==1) fprintf (stdout,"\n number of materials for shrinkage description %ld",numt);
01339 shmat = new shrinkmat [numt];
01340 for (j=0;j<numt;j++){
01341 k=numt+1;
01342 xfscanf (in,"%ld",&k);
01343 if (k>numt || k<1)
01344 print_err("wrong number of materials for shrinkage description",__FILE__,__LINE__,__func__);
01345 shmat[k-1].read (in);
01346 }
01347 break;
01348 }
01349
01350 default:
01351 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
01352 }
01353 }
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367 void mechmat::print (FILE *out)
01368 {
01369 long i, j;
01370
01371 fprintf(out,"\n## materials:\n");
01372 fprintf(out,"\n# number of material types\n");
01373 fprintf (out,"\n%ld\n",nmt);
01374 for (i=0;i<nmt;i++)
01375 {
01376 fprintf (out,"\n\n %d %ld", mtype[i], numtype[i]);
01377 for (j=0;j<numtype[i];j++)
01378 {
01379 fprintf (out,"\n%ld ",j+1);
01380 printmatchar (out, mtype[i], j);
01381 fprintf(out, "\n");
01382 }
01383 }
01384 }
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400 void mechmat::printmatchar (FILE *out, mattype mt, long numinst)
01401 {
01402 switch (mt){
01403 case elisomat:{
01404 eliso[numinst].print (out);
01405 break;
01406 }
01407 case elortomat:{
01408 elorto[numinst].print (out);
01409 break;
01410 }
01411
01412 case scaldamage:{
01413 scdam[numinst].print (out);
01414 break;
01415 }
01416
01417 case nonlocdamgmat:{
01418 nldamg[numinst].print (out);
01419 break;
01420 }
01421
01422 case therisodilat:{
01423 tidilat[numinst].print (out);
01424 break;
01425 }
01426
01427 case relaxationeuro:{
01428 relaxec[numinst].print (out);
01429 break;
01430 }
01431
01432 case shrinkagemat:{
01433 shmat[numinst].print (out);
01434 break;
01435 }
01436
01437 default:
01438 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
01439
01440 }
01441 }
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452 void mechmat::alloceigstrains ()
01453 {
01454 long i,j,nc;
01455
01456 if (eigstrains==NULL){
01457 eigstrains = new double* [tnip];
01458 for (i=0;i<tnip;i++){
01459 nc=ip[i].ncompstr;
01460 eigstrains[i] = new double [nc];
01461 for (j=0;j<nc;j++){
01462 eigstrains[i][j]=0.0;
01463 }
01464 }
01465 }
01466
01467 }
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479 void mechmat::nulleigstrains ()
01480 {
01481 long i,j,nc;
01482
01483 if (eigstrains != NULL){
01484 for (i=0;i<tnip;i++){
01485 nc=ip[i].ncompstr;
01486 for (j=0;j<nc;j++){
01487 eigstrains[i][j]=0.0;
01488 }
01489 }
01490 }
01491 }
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502 void mechmat::alloctempstrains ()
01503 {
01504 long i,j,nc;
01505
01506 if (tempstrains == NULL){
01507 tempstrains = new double* [tnip];
01508 for (i=0;i<tnip;i++){
01509 nc=ip[i].ncompstr;
01510 tempstrains[i] = new double [nc];
01511 for (j=0;j<nc;j++){
01512 tempstrains[i][j]=0.0;
01513 }
01514 }
01515 }
01516
01517 }
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528 void mechmat::nulltempstrains ()
01529 {
01530 long i,j,nc;
01531
01532 if (tempstrains != NULL){
01533 for (i=0;i<tnip;i++){
01534 nc=ip[i].ncompstr;
01535 for (j=0;j<nc;j++){
01536 tempstrains[i][j]=0.0;
01537 }
01538 }
01539 }
01540 }
01541
01542 void mechmat::nullstrains ()
01543 {
01544 long i,j,nc;
01545
01546 for (i=0;i<tnip;i++){
01547 nc=ip[i].ncompstr;
01548 for (j=0;j<nc;j++){
01549 ip[i].strain[j]=0.0;
01550 }
01551 }
01552 }
01553
01554 void mechmat::nullstresses ()
01555 {
01556 long i,j,nc;
01557
01558 for (i=0;i<tnip;i++){
01559 nc=ip[i].ncompstr;
01560 for (j=0;j<nc;j++){
01561 ip[i].stress[j]=0.0;
01562 }
01563 }
01564 }
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574 void mechmat::alloceigstresses ()
01575 {
01576 long i,j,nc;
01577
01578 if (eigstresses == NULL){
01579 eigstresses = new double* [tnip];
01580 for (i=0;i<tnip;i++){
01581 nc=ip[i].ncompstr;
01582 eigstresses[i] = new double [nc];
01583 for (j=0;j<nc;j++){
01584 eigstresses[i][j]=0.0;
01585 }
01586 }
01587 }
01588
01589 }
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602 void mechmat::initmaterialmodels (long lcid)
01603 {
01604 long i,j,nip,ipp,nm;
01605 elemtype te;
01606
01607 for (i=0;i<Mt->ne;i++){
01608 if (Gtm->leso[i]==1){
01609 te = Mt->give_elem_type (i);
01610 nip = Mt->give_tnip (i);
01611
01612 ipp=Mt->elements[i].ipp[0][0];
01613 for(j=0; j<nip; j++){
01614 nm = ip[ipp+j].nm;
01615 initvalues(lcid, ipp+j, 0, 0);
01616 }
01617 }
01618 }
01619 }
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636 void mechmat::initvalues (long lcid, long ipp,long im,long ido)
01637 {
01638 switch (ip[ipp].tm[im])
01639 {
01640 case elisomat:
01641 case elortomat:
01642 case simplas1d:
01643 case jflow:
01644 case jflow2:
01645 case microplaneM4:
01646 case nonlocalmicroM4:
01647 case microsimp:
01648 case microfibro:
01649 case mohrcoul:
01650 case mohrcoulparab:
01651 case boermaterial:
01652 case druckerprager:
01653 case doubledrprager:
01654 case druckerprager2:
01655 break;
01656 case modcamclaymat:
01657 if (Mp->strainstate == 0)
01658 {
01659 compute_ipstrains(lcid);
01660 Mp->strainstate = 1;
01661 }
01662 cclay[ip[ipp].idm[im]].initval(ipp, ido);
01663 break;
01664 case modcamclaycoupmat:
01665 if (Mp->strainstate == 0)
01666 {
01667 compute_ipstrains(lcid);
01668 Mp->strainstate = 1;
01669 }
01670 cclayc[ip[ipp].idm[im]].initval(ipp, ido);
01671 break;
01672 case hypoplastmat:
01673 if (Mp->strainstate == 0)
01674 {
01675 compute_ipstrains(lcid);
01676 Mp->strainstate = 1;
01677 }
01678 hypopl[ip[ipp].idm[im]].initval(ipp, ido);
01679 break;
01680 case shefpl:
01681 case chenplast:
01682 case glasgowmechmat:
01683 case glasgowdamage:
01684 break;
01685 case creep_damage:
01686 crdam[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido);
01687 break;
01688 case shrinkagemat:
01689 shmat[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido);
01690 break;
01691 case time_switchmat:
01692 tswmat[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido);
01693 break;
01694 case effective_stress:
01695 effstr[ip[ipp].idm[im]].initvalues(ipp);
01696 initvalues(lcid,ipp, im+1, ido);
01697 break;
01698 case simvisc:
01699
01700 case lemaitr:
01701 case scaldamage:
01702 case scaldamagecc:
01703 case fixortodamage:
01704 case aci:
01705 case cebfip:
01706 case graphm:
01707 case varelisomat:
01708 case geoelast:
01709 break;
01710 case ortodamage:
01711 ortdam[ip[ipp].idm[im]].initvalues(ipp, ido);
01712 break;
01713 case ortodamagerot:
01714 ortdamrot[ip[ipp].idm[im]].initvalues(ipp, ido);
01715 break;
01716 case anisodamage:
01717 anidam[ip[ipp].idm[im]].initvalues(ipp, ido);
01718 break;
01719 case anisodamagerot:
01720 anidamrot[ip[ipp].idm[im]].initvalues(ipp, ido);
01721 break;
01722 case contmat:
01723 case cebfipcontmat:
01724 break;
01725 case creeprs:
01726 case creepb3:
01727 case creepdpl:{
01728 creep_initmaterialmodel(ipp,im,ido);
01729 }
01730 case creepbaz:
01731 case consolidation:
01732 case winklerpasternak:
01733 case therisodilat:
01734 case relaxationeuro:
01735 case nonlocplastmat:
01736 case nonlocdamgmat:
01737 case damage_plasticity:
01738 case viscoplasticity:
01739 case creepeffym:
01740 case rcmatmodnorm:
01741 break;
01742 default:
01743 print_err("unknown material type is required", __FILE__,__LINE__,__func__);
01744 }
01745 }
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771 void mechmat::storestrain (long lcid,long ipp,vector &eps)
01772 {
01773 long i,j,m,n;
01774 m=lcid*ip[ipp].ncompstr;
01775 n=eps.n;
01776 j=0;
01777 for (i=m;i<m+n;i++){
01778 ip[ipp].strain[i]=eps[j];
01779 j++;
01780 }
01781 }
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797 void mechmat::storestrain (long lcid,long ipp,long fi,vector &eps)
01798 {
01799 long i,j,m,n;
01800 m=lcid*ip[ipp].ncompstr;
01801 n=eps.n;
01802 j=0;
01803 for (i=m+fi;i<m+fi+n;i++){
01804 ip[ipp].strain[i]=eps[j];
01805 j++;
01806 }
01807 }
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825 void mechmat::storestrain (long lcid,long ipp,long fi,long ncomp,vector &eps)
01826 {
01827 long i,j,m;
01828 m=lcid*ip[ipp].ncompstr;
01829 j=fi;
01830 for (i=m+fi;i<m+fi+ncomp;i++){
01831 ip[ipp].strain[i]=eps[j];
01832 j++;
01833 }
01834 }
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852 void mechmat::storestrain (long lcid, long ipp, long fi, long ncomp, double *eps)
01853 {
01854 long i,j,m;
01855 m=lcid*ip[ipp].ncompstr;
01856 j=fi;
01857 for (i=m+fi;i<m+fi+ncomp;i++){
01858 ip[ipp].strain[i]=eps[j];
01859 j++;
01860 }
01861 }
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876 void mechmat::givestrain (long lcid,long ipp,vector &eps)
01877 {
01878 long i,j,m,n;
01879 m=lcid*ip[ipp].ncompstr;
01880 n=eps.n;
01881 j=0;
01882 for (i=m;i<n+m;i++){
01883 eps[j]=ip[ipp].strain[i];
01884 j++;
01885 }
01886 }
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902 void mechmat::givestrain (long lcid,long ipp,long fi,vector &eps)
01903 {
01904 long i,j,m,n;
01905 m=lcid*ip[ipp].ncompstr;
01906 n=eps.n;
01907 j=0;
01908 for (i=m+fi;i<m+fi+n;i++){
01909 eps[j]=ip[ipp].strain[i];
01910 j++;
01911 }
01912 }
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931 void mechmat::givestrain (long lcid,long ipp,long fi,long ncomp,vector &eps)
01932 {
01933 long i,j,m;
01934 m=lcid*ip[ipp].ncompstr;
01935 j=fi;
01936 for (i=m+fi;i<m+fi+ncomp;i++){
01937 eps[j]=ip[ipp].strain[i];
01938 j++;
01939 }
01940 }
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952 void mechmat::cleanstrain ()
01953 {
01954 long i;
01955
01956 for (i=0;i<tnip;i++){
01957
01958 ip[i].clean_strains (1);
01959 }
01960 }
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975 void mechmat::storestress (long lcid,long ipp,vector &sig)
01976 {
01977 long i,j,m,n;
01978 m=lcid*ip[ipp].ncompstr;
01979 n=sig.n;
01980 j=0;
01981 for (i=m;i<m+n;i++){
01982 ip[ipp].stress[i]=sig[j];
01983 j++;
01984 }
01985 }
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001 void mechmat::storestress (long lcid,long ipp,long fi,vector &sig)
02002 {
02003 long i,j,m,n;
02004 m=lcid*ip[ipp].ncompstr;
02005 n=sig.n;
02006 j=0;
02007 for (i=m+fi;i<m+fi+n;i++){
02008 ip[ipp].stress[i]=sig[j];
02009 j++;
02010 }
02011 }
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030 void mechmat::storestress (long lcid,long ipp,long fi,long ncomp,vector &sig)
02031 {
02032 long i,j,m;
02033 m=lcid*ip[ipp].ncompstr;
02034 j=fi;
02035 for (i=m+fi;i<m+fi+ncomp;i++){
02036 ip[ipp].stress[i]=sig[j];
02037 j++;
02038 }
02039 }
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058 void mechmat::storestress (long lcid, long ipp, long fi, long ncomp, double *sig)
02059 {
02060 long i,j,m;
02061 m=lcid*ip[ipp].ncompstr;
02062 j=fi;
02063 for (i=m+fi;i<m+fi+ncomp;i++){
02064 ip[ipp].stress[i]=sig[j];
02065 j++;
02066 }
02067 }
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082 void mechmat::givestress (long lcid,long ipp,vector &sig)
02083 {
02084 long i,j,m,n;
02085 m=lcid*ip[ipp].ncompstr;
02086 n=sig.n;
02087 j=0;
02088 for (i=m;i<m+n;i++){
02089 sig[j]=ip[ipp].stress[i];
02090 j++;
02091 }
02092 }
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106
02107
02108 void mechmat::givestress (long lcid,long ipp,long fi,vector &sig)
02109 {
02110 long i,j,m,n;
02111 m=lcid*ip[ipp].ncompstr;
02112 n=sig.n;
02113 j=0;
02114 for (i=m+fi;i<m+fi+n;i++){
02115 sig[j]=ip[ipp].stress[i];
02116 j++;
02117 }
02118 }
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137 void mechmat::givestress (long lcid,long ipp,long fi,long ncomp,vector &sig)
02138 {
02139 long i,j,m;
02140 m=lcid*ip[ipp].ncompstr;
02141 j=fi;
02142 for (i=m+fi;i<m+fi+ncomp;i++){
02143 sig[j]=ip[ipp].stress[i];
02144 j++;
02145 }
02146 }
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160 void mechmat::storeeigstrain (long ipp,vector &eps)
02161 {
02162 long i,nc;
02163
02164 nc=ip[ipp].ncompstr;
02165 if (est==eigstrain){
02166 for (i=0;i<nc;i++){
02167 eigstrains[ipp][i]=eps[i];
02168 }
02169 }
02170 if (est==tempstrain){
02171 for (i=0;i<nc;i++){
02172 tempstrains[ipp][i]=eps[i];
02173 }
02174 }
02175 }
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191 void mechmat::storeeigstrain (long ipp,long fi,long ncomp,vector &eps)
02192 {
02193 long i,j;
02194
02195 if (est==eigstrain){
02196 j=0;
02197 for (i=fi;i<fi+ncomp;i++){
02198 eigstrains[ipp][i] = eps[j];
02199 j++;
02200 }
02201 }
02202 if (est==tempstrain){
02203 j=0;
02204 for (i=fi;i<fi+ncomp;i++){
02205 tempstrains[ipp][i] = eps[j];
02206 j++;
02207 }
02208 }
02209 }
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221
02222
02223
02224 void mechmat::giveeigstrain (long ipp,vector &eps)
02225 {
02226 long i,n;
02227
02228 n=eps.n;
02229 if (est==eigstrain){
02230 for (i=0;i<n;i++){
02231 eps[i] = eigstrains[ipp][i];
02232 }
02233 }
02234 if (est==tempstrain){
02235 for (i=0;i<n;i++){
02236 eps[i] = tempstrains[ipp][i];
02237 }
02238 }
02239 }
02240
02241
02242
02243
02244
02245
02246
02247
02248
02249
02250
02251
02252
02253
02254
02255
02256 void mechmat::giveeigstrain (long ipp,long fi,long ncomp,vector &eps)
02257 {
02258 long i,j;
02259
02260 if (est==eigstrain){
02261 j=0;
02262 for (i=fi;i<fi+ncomp;i++){
02263 eps[j] = eigstrains[ipp][i];
02264 j++;
02265 }
02266 }
02267 if (est==tempstrain){
02268 j=0;
02269 for (i=fi;i<fi+ncomp;i++){
02270 eps[j] = tempstrains[ipp][i];
02271 j++;
02272 }
02273 }
02274 }
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288 void mechmat::storeeigstress (long ipp,vector &sig)
02289 {
02290 long i,nc;
02291
02292 nc=sig.n;
02293 for (i=0;i<nc;i++){
02294 eigstresses[ipp][i]=sig[i];
02295 }
02296
02297 }
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314 void mechmat::storeeigstress (long ipp,long fi,long ncomp,vector &sig)
02315 {
02316 long i,j;
02317
02318 j=0;
02319 for (i=fi;i<fi+ncomp;i++){
02320 eigstresses[ipp][i] = sig[j];
02321 j++;
02322 }
02323 }
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337 void mechmat::giveeigstress (long ipp,vector &sig)
02338 {
02339 long i,n;
02340
02341 n=sig.n;
02342 for (i=0;i<n;i++){
02343 sig[i] = eigstresses[ipp][i];
02344 }
02345 }
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360 void mechmat::giveeigstress (long ipp,long fi,vector &sig)
02361 {
02362 long i,j,ncomp;
02363 ncomp=sig.n;
02364 j=0;
02365 for (i=fi;i<fi+ncomp;i++){
02366 sig[j] = eigstresses[ipp][i];
02367 j++;
02368 }
02369 }
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385 void mechmat::giveeigstress (long ipp,long fi,long ncomp,vector &sig)
02386 {
02387 long i,j;
02388 j=fi;
02389 for (i=fi;i<fi+ncomp;i++){
02390 sig[j] = eigstresses[ipp][i];
02391 j++;
02392 }
02393 }
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409 void mechmat::storeeqother (long ipp,long fi,long ncomp,double *comp)
02410 {
02411 long i,j;
02412 j=0;
02413 for (i=fi;i<fi+ncomp;i++){
02414 ip[ipp].eqother[i]=comp[j];
02415 j++;
02416 }
02417 }
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433 void mechmat::giveeqother (long ipp,long fi,long ncomp,double *comp)
02434 {
02435 long i,j;
02436 j=0;
02437 for (i=fi;i<fi+ncomp;i++){
02438 comp[j]=ip[ipp].eqother[i];
02439 j++;
02440 }
02441 }
02442
02443
02444
02445
02446
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457 void mechmat::storeother (long ipp,long fi,long ncomp,double *comp)
02458 {
02459 long i,j;
02460 j=0;
02461 for (i=fi;i<fi+ncomp;i++){
02462 ip[ipp].other[i]=comp[j];
02463 j++;
02464 }
02465 }
02466
02467
02468
02469
02470
02471
02472
02473
02474
02475
02476
02477
02478
02479
02480
02481 long mechmat::givencompother (long ipp,long im)
02482 {
02483 long ncompo = 0;
02484 long ncompstr = ip[ipp].ncompstr;
02485
02486 switch (ip[ipp].tm[im])
02487 {
02488 case elisomat:
02489 case elortomat:
02490 case elgmat3d:
02491 case winklerpasternak:
02492 case graphm:
02493 case varelisomat:
02494 case geoelast:
02495 case simplas1d:
02496 case jflow:
02497 case jflow2:
02498 case microplaneM4:
02499 case microsimp:
02500 case microfibro:
02501 case mohrcoul:
02502 case mohrcoulparab:
02503 case boermaterial:
02504 case druckerprager:
02505 case doubledrprager:
02506 case druckerprager2:
02507 case chenplast:
02508 case modcamclaymat:
02509 case modcamclaycoupmat:
02510 case hypoplastmat:
02511 case shefpl:
02512 case scaldamage:
02513 case fixortodamage:
02514 case scaldamagecc:
02515 case ortodamage:
02516 case ortodamagerot:
02517 case anisodamage:
02518 case anisodamagerot:
02519 case creepbaz:
02520 case creepeffym:
02521 case contmat:
02522 case cebfipcontmat:
02523 case consolidation:
02524 case nonlocplastmat:
02525 case nonlocdamgmat:
02526 case nonlocalmicroM4:
02527 case glasgowdamage:
02528 case damage_plasticity:
02529 case therisodilat:
02530 case relaxationeuro:
02531 case simvisc:
02532 case lenjonespot:{
02533 ncompo = givencompeqother (ipp,im);
02534 break;
02535 }
02536 case creeprs:
02537 case creepb3:
02538 case creepdpl:{
02539 ncompo=creep_ncompo(ipp,im);
02540 break;
02541 }
02542 case creep_damage:
02543 ncompo += givencompother (ipp, im+1);
02544 ncompo += givencompother (ipp, im+2);
02545 break;
02546 case shrinkagemat:
02547 ncompo += 1 + 2*ncompstr;
02548 ncompo += givencompother (ipp, im+1);
02549 break;
02550 case effective_stress:
02551 ncompo += givencompother (ipp, im+1);
02552 break;
02553
02554 case lemaitr:
02555 case glasgowmechmat:
02556 case viscoplasticity:
02557 ncompo=0;
02558 break;
02559 case time_switchmat:
02560 ncompo = tswmat[ip[ipp].idm[im]].givencompother (ipp,im);
02561 break;
02562
02563 default:
02564 print_err("unknown material type is required", __FILE__, __LINE__, __func__);
02565 }
02566
02567 return ncompo;
02568 }
02569
02570
02571
02572
02573
02574
02575
02576
02577
02578
02579
02580
02581
02582
02583
02584 long mechmat::givencompeqother (long ipp,long im)
02585 {
02586 long ncompo = 0;
02587 long ncompstr = ip[ipp].ncompstr;
02588 long n;
02589
02590 switch (ip[ipp].tm[im])
02591 {
02592 case elisomat:
02593 break;
02594 case elortomat:
02595 break;
02596 case elgmat3d:
02597 break;
02598
02599
02600
02601 case winklerpasternak:
02602 break;
02603 case graphm:
02604 break;
02605 case varelisomat:
02606 break;
02607 case geoelast:
02608 ncompo = 1;
02609 ncompo += givencompeqother(ipp, im+1);
02610 break;
02611 case simplas1d:
02612 ncompo=3;
02613 break;
02614 case jflow:
02615 ncompo=ncompstr+2;
02616 break;
02617 case jflow2:
02618 ncompo=ncompstr+2;
02619 break;
02620 case microplaneM4:
02621 n=mpM4[ip[ipp].idm[im]].numberOfMicroplanes;
02622 ncompo=3*n+1+ncompstr;
02623 break;
02624 case microsimp:
02625 n=mpSIM[ip[ipp].idm[im]].numberOfMicroplanes;
02626 ncompo=n+1+6;
02627 break;
02628 case microfibro:
02629 n=mpfib[ip[ipp].idm[im]].numberOfMicroplanes;
02630 ncompo=5*n+1+6;
02631 break;
02632 case mohrcoul:
02633 ncompo=ncompstr+1;
02634 break;
02635 case mohrcoulparab:
02636 ncompo=ncompstr+1;
02637 break;
02638 case boermaterial:
02639 ncompo=ncompstr+1;
02640 break;
02641 case doubledrprager:
02642 ncompo=ncompstr+4;
02643 break;
02644 case druckerprager:
02645 ncompo=ncompstr+4;
02646 break;
02647 case druckerprager2:
02648 ncompo=ncompstr+4;
02649 break;
02650 case chenplast:
02651 ncompo=ncompstr+3;
02652 break;
02653 case modcamclaymat:
02654 ncompo=ncompstr+2+3+4;
02655 break;
02656 case modcamclaycoupmat:
02657 ncompo=ncompstr+3+2+4;
02658 case hypoplastmat:
02659 ncompo=2*ncompstr+hypopl[ip[ipp].idm[im]].nstatv+1;
02660 break;
02661 case shefpl:
02662 ncompo=ncompstr+2;
02663 break;
02664 case glasgowmechmat:
02665 ncompo=ncompstr*4+5;
02666 break;
02667 case glasgowdamage:
02668 ncompo=2;
02669 break;
02670 case creep_damage:
02671 ncompo += givencompeqother (ipp, im+1);
02672 ncompo += givencompeqother (ipp, im+2);
02673 break;
02674 case shrinkagemat:
02675 ncompo += 1 + 2*ncompstr;
02676 ncompo += givencompeqother (ipp, im+1);
02677 break;
02678 case time_switchmat:
02679 ncompo = tswmat[ip[ipp].idm[im]].givencompeqother (ipp,im);
02680 break;
02681 case effective_stress:
02682 ncompo += givencompeqother (ipp, im+1);
02683 break;
02684 case simvisc:
02685 ncompo=3*ncompstr;
02686 break;
02687
02688
02689
02690 case lemaitr:
02691 ncompo=3*ncompstr+2;
02692 break;
02693 case scaldamage:
02694 ncompo=3;
02695 break;
02696 case fixortodamage:
02697 ncompo=3;
02698 break;
02699 case scaldamagecc:
02700 ncompo=2;
02701 break;
02702 case ortodamage:
02703 ncompo=19;
02704 break;
02705 case ortodamagerot:
02706 ncompo=23;
02707 break;
02708 case anisodamage:
02709
02710 ncompo=26;
02711 break;
02712 case anisodamagerot:
02713
02714 ncompo=22;
02715 break;
02716 case cebfipcontmat:
02717 case contmat:
02718 ncompo = 2;
02719 break;
02720 case creeprs:
02721 case creepb3:
02722 case creepdpl:{
02723 ncompo=creep_ncompo(ipp,im);
02724 break;
02725 }
02726 case creepbaz:
02727
02728 n=crbaz[ip[ipp].idm[im]].numberOfCreepb();
02729 ncompo=ncompstr*(n+3)+3+1;
02730 break;
02731 case creepeffym:
02732 ncompo=ncompstr;
02733 break;
02734 case consolidation:
02735 n=csol[ip[ipp].idm[im]].numberOfConsol();
02736 if(ncompstr==6) ncompo=ncompstr+2*n+1;
02737 if(ncompstr==5) ncompo=ncompstr+3*n+1;
02738 if(ncompstr==3) ncompo=ncompstr+3*n+1;
02739
02740
02741 break;
02742
02743 case nonlocplastmat:
02744 ncompo = givencompeqother (ipp, im+1);
02745 break;
02746 case nonlocdamgmat:
02747 ncompo = givencompeqother (ipp, im+1);
02748 break;
02749 case nonlocalmicroM4:
02750 n=mpM4[ip[ipp].idm[im]].numberOfMicroplanes;
02751 ncompo=3*n+1+ncompstr+ncompstr;
02752 break;
02753 case damage_plasticity:
02754 ncompo = givencompeqother (ipp, im+1);
02755 ncompo += givencompeqother (ipp, im+2);
02756 break;
02757 case viscoplasticity:
02758 ncompo = givencompeqother (ipp, im+1);
02759 ncompo += givencompeqother (ipp, im+2);
02760 break;
02761 case therisodilat:{
02762 break;
02763 }
02764 case relaxationeuro:{
02765 break;
02766 }
02767 case lenjonespot:{
02768 ncompo = 0;
02769 break;
02770 }
02771 default:{
02772 print_err("unknown material type is required", __FILE__, __LINE__, __func__);
02773 }
02774 }
02775
02776 if (((ip[ipp].hmt & 1) && (im < ip[ipp].nm-1)) && (im == 0))
02777
02778
02779
02780 {
02781 switch (ip[ipp].tm[ip[ipp].nm-1])
02782 {
02783 case therisodilat:
02784 break;
02785 default:{
02786 print_err("unknown thermal material type is required", __FILE__, __LINE__, __func__);
02787 }
02788 }
02789 }
02790 return ncompo;
02791 }
02792
02793
02794
02795
02796
02797
02798
02799
02800
02801
02802
02803
02804
02805
02806
02807 void mechmat::storenonloc (long ipp,long fi,long ncomp,double *comp)
02808 {
02809 long i,j;
02810 j=0;
02811 for (i=fi;i<fi+ncomp;i++){
02812 ip[ipp].nonloc[i]=comp[j];
02813 j++;
02814 }
02815 }
02816
02817
02818
02819
02820
02821
02822
02823
02824
02825
02826
02827
02828
02829
02830
02831 void mechmat::givenonloc (long ipp,long fi,long ncomp,double *comp)
02832 {
02833 long i,j;
02834 j=0;
02835 for (i=fi;i<fi+ncomp;i++){
02836 comp[j]=ip[ipp].nonloc[i];
02837 j++;
02838 }
02839 }
02840
02841
02842
02843
02844
02845
02846
02847
02848
02849
02850
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869 void mechmat::givequantity (integratedquant iq,long lcid,long ipp,long fi,vector &ipv)
02870 {
02871 switch (iq){
02872 case locstress:{
02873
02874 givestress (lcid,ipp,fi,ipv);
02875 break;
02876 }
02877 case nonlocstress:{
02878
02879 givestress (lcid,ipp,fi,ipv);
02880 break;
02881 }
02882 case stressincr:{
02883
02884 givestressincr (lcid,ipp,0,0,fi,ipv);
02885 break;
02886 }
02887 case eigstress:{
02888
02889 giveeigstress (ipp,fi,ipv);
02890 break;
02891 }
02892 default:{
02893 print_err("unknown type of quantity is required",__FILE__,__LINE__,__func__);
02894 }
02895 }
02896 }
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912 double mechmat::give_actual_ym(long ipp, long im, long ido)
02913 {
02914 long ncompo;
02915 double e = 0.0;
02916
02917 switch (ip[ipp].tm[im])
02918 {
02919 case elisomat:
02920 e = eliso[ip[ipp].idm[im]].e;
02921 break;
02922 case modcamclaymat:
02923 e = cclay[ip[ipp].idm[im]].give_actual_ym(ipp, im, ido);
02924 break;
02925 case winklerpasternak:
02926 case simplas1d:
02927 case jflow:
02928 case jflow2:
02929 case microplaneM4:
02930 case microsimp:
02931 case microfibro:
02932 case mohrcoul:
02933 case mohrcoulparab:
02934 case boermaterial:
02935 case druckerprager:
02936 case doubledrprager:
02937 case druckerprager2:
02938 case chenplast:
02939 case shefpl:
02940 case glasgowmechmat:
02941 case glasgowdamage:
02942 case simvisc:
02943
02944 case lemaitr:
02945 case scaldamage:
02946 case fixortodamage:
02947 case scaldamagecc:
02948 case ortodamage:
02949 case ortodamagerot:
02950 case anisodamage:
02951 case anisodamagerot:
02952 case creepdpl:
02953 case creepbaz:
02954 case creepeffym:
02955 case consolidation:
02956 ncompo = givencompeqother(ipp, im);
02957 e = give_actual_ym(ipp, im+1, ido+ncompo);
02958 break;
02959 case damage_plasticity:
02960 ncompo = givencompeqother(ipp, im+1);
02961 e= give_actual_ym(ipp, im+2, ido+ncompo);
02962 break;
02963 case viscoplasticity:
02964 ncompo = givencompeqother(ipp, im+1);
02965 e = give_actual_ym(ipp, im+2, ido+ncompo);
02966 break;
02967 case nonlocplastmat:
02968 case nonlocdamgmat:
02969 case nonlocalmicroM4:
02970 case creep_damage:
02971 e = give_actual_ym(ipp, im+1, ido);
02972 break;
02973 case shrinkagemat:
02974 ncompo = givencompeqother(ipp, im);
02975 e = give_actual_ym (ipp, im+1, ido+ncompo);
02976 break;
02977 case creepb3:
02978 case creeprs:
02979 e = creep_give_actual_ym(ipp, im, ido);
02980 break;
02981 case time_switchmat:
02982 e = tswmat[ip[ipp].idm[im]].give_actual_ym(ipp, im, ido);
02983 break;
02984 case effective_stress:{
02985 ncompo = givencompeqother(ipp, im);
02986 e = give_actual_ym(ipp, im+1, ido+ncompo);
02987 break;
02988 }
02989 default:
02990 print_err("unknown elastic material type is required",__FILE__,__LINE__,__func__);
02991 }
02992
02993 return e;
02994 }
02995
02996
02997
02998
02999
03000
03001
03002
03003
03004
03005
03006
03007
03008
03009
03010 double mechmat::give_initial_ym(long ipp, long im, long ido)
03011 {
03012 long ncompo;
03013 double e = 0.0;
03014
03015 switch (ip[ipp].tm[im])
03016 {
03017 case elisomat:
03018 e = eliso[ip[ipp].idm[im]].e;
03019 break;
03020 case varelisomat:
03021 e = veliso[ip[ipp].idm[im]].e;
03022 break;
03023 case winklerpasternak:
03024 case simplas1d:
03025 case jflow:
03026 case jflow2:
03027 case microplaneM4:
03028 case microsimp:
03029 case microfibro:
03030 case mohrcoul:
03031 case mohrcoulparab:
03032 case boermaterial:
03033 case druckerprager:
03034 case doubledrprager:
03035 case druckerprager2:
03036 case chenplast:
03037 case modcamclaymat:
03038 case shefpl:
03039 case glasgowmechmat:
03040 case glasgowdamage:
03041 case simvisc:
03042
03043 case lemaitr:
03044 case scaldamage:
03045 case fixortodamage:
03046 case scaldamagecc:
03047 case ortodamage:
03048 case ortodamagerot:
03049 case anisodamage:
03050 case anisodamagerot:
03051 case creepdpl:
03052 case creepbaz:
03053 case creepeffym:
03054 case consolidation:
03055 ncompo = givencompeqother(ipp, im);
03056 e = give_initial_ym(ipp, im+1, ido+ncompo);
03057 break;
03058 case damage_plasticity:
03059 ncompo = givencompeqother(ipp, im+1);
03060 e= give_initial_ym(ipp, im+2, ido+ncompo);
03061 break;
03062 case viscoplasticity:
03063 ncompo = givencompeqother(ipp, im+1);
03064 e = give_initial_ym(ipp, im+2, ido+ncompo);
03065 break;
03066 case nonlocplastmat:
03067 case nonlocdamgmat:
03068 case nonlocalmicroM4:
03069 case creep_damage:
03070 case effective_stress:
03071 e = give_initial_ym(ipp, im+1, ido);
03072 break;
03073 case shrinkagemat:
03074 ncompo = givencompeqother(ipp, im);
03075 e = give_initial_ym (ipp, im+1, ido+ncompo);
03076 break;
03077 case creepb3:
03078 case creeprs:
03079 e = creep_compute_inital_ym(ipp, im, ido);
03080 break;
03081 case time_switchmat:
03082 e = tswmat[ip[ipp].idm[im]].give_initial_ym(ipp, im, ido);
03083 break;
03084 default:
03085 print_err("unknown elastic material type is required",__FILE__,__LINE__,__func__);
03086 }
03087 return e;
03088 }
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103
03104 double mechmat::give_actual_nu(long ipp, long im, long ido)
03105 {
03106 long ncompo;
03107 double nu = 0.0;
03108
03109 switch (ip[ipp].tm[im])
03110 {
03111 case elisomat:
03112 nu = eliso[ip[ipp].idm[im]].nu;
03113 break;
03114 case modcamclaymat:
03115 nu = cclay[ip[ipp].idm[im]].give_actual_nu(ipp, im, ido);
03116 break;
03117 case winklerpasternak:
03118 case simplas1d:
03119 case jflow:
03120 case jflow2:
03121 case microplaneM4:
03122 case microsimp:
03123 case microfibro:
03124 case mohrcoul:
03125 case mohrcoulparab:
03126 case boermaterial:
03127 case druckerprager:
03128 case doubledrprager:
03129 case druckerprager2:
03130 case chenplast:
03131 case shefpl:
03132 case glasgowmechmat:
03133 case glasgowdamage:
03134 case simvisc:
03135
03136 case lemaitr:
03137 case scaldamage:
03138 case fixortodamage:
03139 case scaldamagecc:
03140 case ortodamage:
03141 case ortodamagerot:
03142 case anisodamage:
03143 case anisodamagerot:
03144 case creepdpl:
03145 case creepbaz:
03146 case creepb3:
03147 case creeprs:
03148 case creepeffym:
03149 case consolidation:
03150 ncompo = givencompeqother(ipp, im);
03151 nu = give_actual_nu(ipp, im+1, ido+ncompo);
03152 break;
03153 case damage_plasticity:
03154 ncompo = givencompeqother(ipp, im+1);
03155 nu = give_actual_nu(ipp, im+2, ido+ncompo);
03156 break;
03157 case viscoplasticity:
03158 ncompo = givencompeqother(ipp, im+1);
03159 nu = give_actual_nu(ipp, im+2, ido+ncompo);
03160 break;
03161 case creep_damage:
03162 ncompo = givencompeqother(ipp, im+1);
03163 nu = give_actual_nu(ipp, im+2, ido+ncompo);
03164 break;
03165 case shrinkagemat:
03166 ncompo = givencompeqother(ipp, im);
03167 nu = give_actual_nu(ipp, im+1, ido+ncompo);
03168 break;
03169 case nonlocplastmat:
03170 case nonlocdamgmat:
03171 case nonlocalmicroM4:
03172 case effective_stress:
03173 nu = give_actual_nu(ipp, im+1, ido);
03174 break;
03175 case time_switchmat:
03176 nu = tswmat[ip[ipp].idm[im]].give_actual_nu(ipp, im, ido);
03177 break;
03178 default:
03179 print_err("unknown elastic material type is required",__FILE__,__LINE__,__func__);
03180 }
03181 return nu;
03182 }
03183
03184
03185
03186
03187
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
03198 double mechmat::give_initial_nu(long ipp, long im, long ido)
03199 {
03200 long ncompo;
03201 double nu = 0.0;
03202
03203 switch (ip[ipp].tm[im])
03204 {
03205 case elisomat:
03206 nu = eliso[ip[ipp].idm[im]].nu;
03207 break;
03208 case varelisomat:
03209 nu = veliso[ip[ipp].idm[im]].nu;
03210 break;
03211 case winklerpasternak:
03212 case simplas1d:
03213 case jflow:
03214 case jflow2:
03215 case microplaneM4:
03216 case microsimp:
03217 case microfibro:
03218 case mohrcoul:
03219 case mohrcoulparab:
03220 case boermaterial:
03221 case druckerprager:
03222 case doubledrprager:
03223 case druckerprager2:
03224 case chenplast:
03225 case modcamclaymat:
03226 case shefpl:
03227 case glasgowmechmat:
03228 case glasgowdamage:
03229 case simvisc:
03230
03231 case lemaitr:
03232 case scaldamage:
03233 case fixortodamage:
03234 case scaldamagecc:
03235 case ortodamage:
03236 case ortodamagerot:
03237 case anisodamage:
03238 case anisodamagerot:
03239 case creepdpl:
03240 case creepbaz:
03241 case creepb3:
03242 case creeprs:
03243 case creepeffym:
03244 case consolidation:
03245 ncompo = givencompeqother(ipp, im);
03246 nu = give_actual_nu(ipp, im+1, ido+ncompo);
03247 break;
03248 case damage_plasticity:
03249 ncompo = givencompeqother(ipp, im+1);
03250 nu = give_actual_nu(ipp, im+2, ido+ncompo);
03251 break;
03252 case viscoplasticity:
03253 ncompo = givencompeqother(ipp, im+1);
03254 nu = give_actual_nu(ipp, im+2, ido+ncompo);
03255 break;
03256 case creep_damage:
03257 ncompo = givencompeqother(ipp, im+1);
03258 nu = give_actual_nu(ipp, im+2, ido+ncompo);
03259 break;
03260 case shrinkagemat:
03261 ncompo = givencompeqother(ipp, im);
03262 nu = give_actual_nu(ipp, im+1, ido+ncompo);
03263 break;
03264 case nonlocplastmat:
03265 case nonlocdamgmat:
03266 case nonlocalmicroM4:
03267 case effective_stress:
03268 nu = give_actual_nu(ipp, im+1, ido);
03269 break;
03270 case time_switchmat:
03271 nu = tswmat[ip[ipp].idm[im]].give_actual_nu(ipp, im, ido);
03272 break;
03273 default:
03274 print_err("unknown elastic material type is required",__FILE__,__LINE__,__func__);
03275 }
03276 return nu;
03277 }
03278
03279
03280
03281
03282
03283
03284
03285
03286
03287
03288
03289
03290
03291
03292
03293 double mechmat::give_actual_ft(long ipp, long im, long ido)
03294 {
03295 double ft = 0.0;
03296
03297 switch (ip[ipp].tm[im])
03298 {
03299 case scaldamage:
03300 ft = scdam[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
03301 break;
03302 case ortodamage:
03303 ft = ortdam[ip[ipp].idm[im]].give_actual_ft(ipp);
03304 break;
03305 case ortodamagerot:
03306 ft = ortdamrot[ip[ipp].idm[im]].give_actual_ft(ipp);
03307 break;
03308 case anisodamage:
03309 ft = anidam[ip[ipp].idm[im]].give_actual_ft(ipp);
03310 break;
03311 case anisodamagerot:
03312 ft = anidamrot[ip[ipp].idm[im]].give_actual_ft(ipp);
03313 break;
03314 case nonlocdamgmat:
03315 ft = give_actual_ft(ipp,im+1,ido);
03316 break;
03317 case damage_plasticity:
03318 ft = dampl[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
03319 break;
03320 case creep_damage:
03321 ft = crdam[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
03322 break;
03323 case shrinkagemat:
03324 ft = shmat[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
03325 break;
03326 case effective_stress:
03327 ft = give_actual_ft(ipp,im+1,ido);
03328 break;
03329 case creepb3:
03330 case creeprs:
03331 ft = creep_give_actual_ft(ipp, im, ido);
03332 break;
03333 case creepeffym:
03334 ft = creffym[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
03335 break;
03336 case time_switchmat:
03337 ft = tswmat[ip[ipp].idm[im]].give_actual_ft(ipp, im, ido);
03338 break;
03339 default:
03340 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
03341 }
03342 return ft;
03343 }
03344
03345
03346
03347
03348
03349
03350
03351
03352
03353
03354
03355
03356
03357
03358
03359 double mechmat::give_actual_fc(long ipp, long im, long ido)
03360 {
03361 double fc = 0.0;
03362
03363 switch (ip[ipp].tm[im])
03364 {
03365 case ortodamage:
03366 fc = ortdam[ip[ipp].idm[im]].give_actual_fc(ipp);
03367 break;
03368 case ortodamagerot:
03369 fc = ortdamrot[ip[ipp].idm[im]].give_actual_fc(ipp);
03370 break;
03371 case anisodamage:
03372 fc = anidam[ip[ipp].idm[im]].give_actual_fc(ipp);
03373 break;
03374 case anisodamagerot:
03375 fc = anidamrot[ip[ipp].idm[im]].give_actual_fc(ipp);
03376 break;
03377 case nonlocdamgmat:
03378 fc = give_actual_fc(ipp,im+1,ido);
03379 break;
03380 case creep_damage:
03381 fc = crdam[ip[ipp].idm[im]].give_actual_fc(ipp, im, ido);
03382 break;
03383 case shrinkagemat:
03384 fc = shmat[ip[ipp].idm[im]].give_actual_fc(ipp, im, ido);
03385 break;
03386 default:
03387 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
03388 }
03389 return fc;
03390 }
03391
03392
03393
03394
03395
03396
03397
03398
03399
03400
03401
03402
03403
03404
03405 double mechmat::givemechq(nontransquant qt, long ipp)
03406 {
03407 double ret=0.0;
03408
03409 switch (qt){
03410 case precons_press:{
03411
03412 ret = give_preconspress(ipp);
03413 break;
03414 }
03415 case mean_stress_eff:{
03416 vector sig(Mm->ip[ipp].ncompstr);
03417 matrix sigt(3,3);
03418 givestress(0,ipp,sig);
03419 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
03420
03421 ret = first_invar(sigt)/3.0;
03422 if (Mm->givestatusnmq(pore_pressure))
03423 {
03424
03425 ret -= Mm->givenonmechq(pore_pressure, ipp);
03426 }
03427 break;
03428 }
03429 case virgin_porosity:{
03430
03431 ret = give_virgporosity(ipp);
03432 break;
03433 }
03434 case init_porosity:{
03435
03436 ret = give_iniporosity(ipp);
03437 break;
03438 }
03439 case scal_iso_damage:{
03440
03441 ret = give_dampar (ipp);
03442 break;
03443 }
03444 case proc_zone_length:{
03445
03446 ret = give_proczonelength (ipp);
03447 break;
03448 }
03449 case crack_width:{
03450
03451 ret = give_crackwidth (ipp);
03452 break;
03453 }
03454 default:{
03455 print_err("unknown type of quantity is required",__FILE__,__LINE__,__func__);
03456 }
03457 }
03458
03459 return ret;
03460 }
03461
03462
03463
03464
03465
03466
03467
03468
03469
03470
03471
03472
03473
03474
03475
03476 double mechmat::give_preconspress(long ipp, long im, long ido)
03477 {
03478 long i;
03479 double pc;
03480
03481 i = ip[ipp].idm[im];
03482
03483 switch (ip[ipp].tm[im]){
03484 case modcamclaymat:{
03485 pc = cclay[i].give_preconspress(ipp,ido);
03486 break;
03487 }
03488 case modcamclaycoupmat:{
03489 pc = cclayc[i].give_preconspress(ipp,ido);
03490 break;
03491 }
03492 case effective_stress:
03493 pc = give_preconspress(ipp, im+1, ido);
03494 break;
03495 case damage_plasticity:
03496 pc = give_preconspress(ipp, im+1, ido);
03497 break;
03498 default:{
03499 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
03500 }
03501 }
03502
03503 return pc;
03504 }
03505
03506
03507
03508
03509
03510
03511
03512
03513
03514
03515
03516
03517
03518
03519
03520 double mechmat::give_virgporosity(long ipp, long im, long ido)
03521 {
03522 long i;
03523 double e_lambda1;
03524
03525 i = ip[ipp].idm[im];
03526
03527 switch (ip[ipp].tm[im]){
03528 case modcamclaymat:{
03529 e_lambda1 = cclay[i].give_virgporosity(ipp,ido);
03530 break;
03531 }
03532 case modcamclaycoupmat:{
03533 e_lambda1 = cclayc[i].give_virgporosity(ipp,ido);
03534 break;
03535 }
03536 case hypoplastmat:{
03537 e_lambda1 = hypopl[i].give_virgporosity(ipp,ido);
03538 break;
03539 }
03540 case effective_stress:
03541 e_lambda1 = give_virgporosity(ipp, im+1, ido);
03542 break;
03543 case damage_plasticity:
03544 e_lambda1 = give_virgporosity(ipp, im+1, ido);
03545 break;
03546 default:{
03547 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
03548 }
03549 }
03550
03551 return e_lambda1;
03552 }
03553
03554
03555
03556
03557
03558
03559
03560
03561
03562
03563
03564
03565
03566
03567
03568 double mechmat::give_iniporosity(long ipp, long im, long ido)
03569 {
03570 long i;
03571 double e_ini;
03572
03573 i = ip[ipp].idm[im];
03574
03575 switch (ip[ipp].tm[im]){
03576 case modcamclaymat:{
03577 e_ini = cclay[i].give_iniporosity(ipp,ido);
03578 break;
03579 }
03580 case modcamclaycoupmat:{
03581 e_ini = cclayc[i].give_iniporosity(ipp,ido);
03582 break;
03583 }
03584 case hypoplastmat:{
03585 e_ini = hypopl[i].give_iniporosity(ipp,ido);
03586 break;
03587 }
03588 case effective_stress:
03589 e_ini = give_iniporosity(ipp, im+1, ido);
03590 break;
03591 case damage_plasticity:
03592 e_ini = give_iniporosity(ipp, im+1, ido);
03593 break;
03594 default:{
03595 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
03596 }
03597 }
03598
03599 return e_ini;
03600 }
03601
03602
03603
03604
03605
03606
03607
03608
03609
03610
03611
03612
03613
03614
03615 double mechmat::give_consparam (long ipp, long im, long ido)
03616 {
03617 long i;
03618 double cp;
03619
03620 i = ip[ipp].idm[im];
03621
03622 switch (ip[ipp].tm[im]){
03623 case simplas1d:{
03624 cp = spl1d[i].give_consparam (ipp,ido);
03625 break;
03626 }
03627 case jflow:{
03628 cp = j2f[i].give_consparam (ipp,ido);
03629 break;
03630 }
03631 case mohrcoul:{
03632 cp = mcoul[i].give_consparam (ipp,ido);
03633 break;
03634 }
03635 case mohrcoulparab:{
03636 cp = mcpar[i].give_consparam (ipp,ido);
03637 break;
03638 }
03639 case boermaterial:{
03640 cp = boerm[i].give_consparam (ipp,ido);
03641 break;
03642 }
03643 case druckerprager:{
03644 cp = drprm[i].give_consparam (ipp,ido);
03645 break;
03646 }
03647 case doubledrprager:{
03648
03649 break;
03650 }
03651 case druckerprager2:{
03652 cp = drprm2[i].give_consparam (ipp,ido);
03653 break;
03654 }
03655 case modcamclaymat:{
03656 cp = cclay[i].give_consparam (ipp,ido);
03657 break;
03658 }
03659 case modcamclaycoupmat:{
03660 cp = cclayc[i].give_consparam (ipp,ido);
03661 break;
03662 }
03663 case hissplasmat:{
03664 cp = hisspl[i].give_consparam (ipp,ido);
03665 break;
03666 }
03667 case damage_plasticity:
03668 cp = give_consparam(ipp, im+1, ido);
03669 break;
03670 case effective_stress:
03671 cp = give_consparam(ipp, im+1, ido);
03672 break;
03673 default:{
03674 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
03675 }
03676 }
03677
03678 return cp;
03679 }
03680
03681
03682
03683
03684
03685
03686
03687
03688
03689
03690
03691
03692
03693
03694 double mechmat::give_dampar(long ipp, long im, long ido)
03695 {
03696 long tncompo;
03697 long damncompo;
03698 double dam = 0.0;
03699 switch (ip[ipp].tm[im])
03700 {
03701 case elisomat:
03702 dam = 0.0;
03703 break;
03704 case elortomat:
03705 dam = 0.0;
03706 break;
03707 case scaldamage:{
03708 dam = scdam[ip[ipp].idm[im]].givedamage(ipp, ido);
03709 break;
03710 }
03711 case nonlocdamgmat:{
03712 dam = give_dampar(ipp, im+1, ido);
03713 break;
03714 }
03715 case shrinkagemat:{
03716 tncompo = givencompeqother(ipp,im);
03717 damncompo = givencompeqother(ipp,im+1);
03718 dam = give_dampar(ipp, im+1, ido+tncompo-damncompo);
03719 break;
03720 }
03721 case effective_stress:{
03722 dam = give_dampar (ipp, im+1, ido);
03723 break;
03724 }
03725 default:{
03726 print_err("unknown material model is required",__FILE__,__LINE__,__func__);
03727 }
03728 }
03729 return dam;
03730 }
03731
03732
03733
03734
03735
03736
03737
03738
03739
03740
03741
03742
03743 double mechmat::give_proczonelength (long ipp, long im, long ido)
03744 {
03745 long tncompo;
03746 long damncompo;
03747 double l = 0.0;
03748 switch (ip[ipp].tm[im])
03749 {
03750 case elisomat:
03751 l = 0.0;
03752 break;
03753 case elortomat:
03754 l = 0.0;
03755 break;
03756 case scaldamage:{
03757 l = scdam[ip[ipp].idm[im]].give_proczonelength(ipp, ido);
03758 break;
03759 }
03760 case effective_stress:{
03761 l = give_proczonelength (ipp, im+1, ido);
03762 break;
03763 }
03764 case nonlocdamgmat:
03765 l = nldamg[ip[ipp].idm[im]].give_proczonelength (ipp, im, ido);
03766 break;
03767 case shrinkagemat:
03768 tncompo = givencompeqother(ipp,im);
03769 damncompo = givencompeqother(ipp,im+1);
03770 l = give_proczonelength (ipp, im+1, ido+tncompo-damncompo);
03771 break;
03772 default:{
03773 print_err("unknown material model is required",__FILE__,__LINE__,__func__);
03774 }
03775 }
03776 return l;
03777 }
03778
03779
03780
03781
03782
03783
03784
03785
03786
03787
03788
03789
03790
03791 double mechmat::give_crackwidth (long ipp, long im, long ido)
03792 {
03793 long tncompo;
03794 long damncompo;
03795 double dam = 0.0;
03796 switch (ip[ipp].tm[im])
03797 {
03798 case elisomat:
03799 dam = 0.0;
03800 break;
03801 case elortomat:
03802 dam = 0.0;
03803 break;
03804 case scaldamage:{
03805 dam = scdam[ip[ipp].idm[im]].give_crackwidth (ipp, ido);
03806 break;
03807 }
03808 case effective_stress:{
03809 dam = give_crackwidth (ipp, im+1, ido);
03810 break;
03811 }
03812 case nonlocdamgmat:
03813 dam = give_crackwidth (ipp, im+1, ido);
03814 break;
03815 case shrinkagemat:
03816 tncompo = givencompeqother(ipp,im);
03817 damncompo = givencompeqother(ipp,im+1);
03818 dam = give_crackwidth (ipp, im+1, ido+tncompo-damncompo);
03819 break;
03820 default:{
03821 print_err("unknown material model is required",__FILE__,__LINE__,__func__);
03822 }
03823 }
03824 return dam;
03825 }
03826
03827
03828
03829
03830
03831
03832
03833
03834
03835
03836
03837
03838
03839
03840
03841 void mechmat::give_kappa(long ipp,long im,long ido, vector &kappa)
03842 {
03843 switch (ip[ipp].tm[im])
03844 {
03845 case scaldamage:
03846 kappa[0] = ip[ipp].other[ido];
03847 break;
03848 case fixortodamage:
03849 kappa[0] = ip[ipp].other[ido];
03850 break;
03851 case scaldamagecc:
03852 kappa[0] = ip[ipp].other[ido];
03853
03854 break;
03855 default:
03856 print_err ("unknown type of damage model is required",__FILE__,__LINE__,__func__);
03857 }
03858 return;
03859 }
03860
03861
03862
03863
03864
03865
03866
03867
03868
03869
03870
03871
03872
03873
03874
03875 void mechmat::giveirrstrains(long ipp, long im, long ido, vector &epsirr)
03876 {
03877 fillv(0.0, epsirr);
03878 switch (ip[ipp].tm[im])
03879 {
03880 case simplas1d:{
03881 spl1d[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03882 break;
03883 }
03884 case jflow:{
03885 j2f[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03886 break;
03887 }
03888 case jflow2:{
03889 j2f2[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03890 break;
03891 }
03892 case mohrcoul:{
03893 mcoul[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03894 break;
03895 }
03896 case mohrcoulparab:{
03897 mcpar[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03898 break;
03899 }
03900 case boermaterial:{
03901 boerm[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03902 break;
03903 }
03904 case druckerprager:{
03905 drprm[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03906 break;
03907 }
03908 case doubledrprager:{
03909
03910 break;
03911 }
03912 case druckerprager2:{
03913 drprm2[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03914 break;
03915 }
03916 case chenplast:{
03917
03918 break;
03919 }
03920 case modcamclaymat:{
03921 cclay[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03922 break;
03923 }
03924 case modcamclaycoupmat:{
03925 cclayc[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03926 break;
03927 }
03928 case hissplasmat:{
03929 hisspl[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03930 break;
03931 }
03932 case effective_stress:{
03933 giveirrstrains(ipp, im+1, ido, epsirr);
03934 break;
03935 }
03936 case creepbaz:{
03937 crbaz[ip[ipp].idm[im]].giveirrstrains(ipp, ido, epsirr);
03938 break;
03939 }
03940 case creepeffym:{
03941 giveirrstrains(ipp, im+1, ido, epsirr);
03942 break;
03943 }
03944 case creepb3:
03945 case creeprs:
03946 case creepdpl:{
03947 creep_giveirrstrains(ipp, im, ido, epsirr);
03948 break;
03949 }
03950 case shrinkagemat:
03951 shmat[ip[ipp].idm[im]].giveirrstrains(ipp, im, ido, epsirr);
03952 default:
03953 print_err("unknown material model is required",__FILE__,__LINE__,__func__);
03954 }
03955 }
03956
03957
03958
03959
03960
03961
03962
03963
03964
03965
03966
03967
03968
03969
03970
03971
03972
03973
03974
03975
03976
03977
03978
03979
03980
03981
03982
03983
03984
03985
03986 long mechmat::search_reqnmq(nonmechquant* &rnmq)
03987 {
03988 long i, j, n;
03989 long anmq[tnknmq];
03990 memset(anmq, 0, sizeof(*anmq)*tnknmq);
03991
03992
03993 for (i=0; i<tnip; i++)
03994 {
03995 for (j=0; j<ip[i].nm; j++)
03996 give_reqnmq(i, j, anmq);
03997 }
03998
03999
04000 n = 0;
04001 for (i=0; i<tnknmq; i++)
04002 {
04003 if (anmq[i] == 1)
04004 n++;
04005 }
04006
04007
04008 if (n == 0)
04009 {
04010 rnmq = NULL;
04011 return n;
04012 }
04013 else
04014 {
04015
04016 rnmq = new nonmechquant[n];
04017 j = 0;
04018 for (i=0; i<tnknmq; i++)
04019 {
04020 if (anmq[i] == 1){
04021 rnmq[j] = nonmechquant(i+1);
04022 j++;
04023 }
04024 }
04025 }
04026
04027 return n;
04028 }
04029
04030
04031
04032
04033
04034
04035
04036
04037
04038
04039
04040
04041
04042
04043
04044
04045
04046
04047
04048 void mechmat::give_reqnmq(long ipp, long im, long *anmq)
04049 {
04050 switch(ip[ipp].tm[im])
04051 {
04052 case elisomat:
04053 case elortomat:
04054 case elgmat3d:
04055 case winklerpasternak:
04056 case graphm:
04057 case varelisomat:
04058 case geoelast:
04059 case simplas1d:
04060 case jflow:
04061 case jflow2:
04062 case microplaneM4:
04063 case microsimp:
04064 case microfibro:
04065 case mohrcoul:
04066 case mohrcoulparab:
04067 case boermaterial:
04068 case druckerprager:
04069 case doubledrprager:
04070 case druckerprager2:
04071 case chenplast:
04072 case modcamclaymat:
04073 case shefpl:
04074 case scaldamage:
04075 case fixortodamage:
04076 case scaldamagecc:
04077 case ortodamage:
04078 case ortodamagerot:
04079 case anisodamage:
04080 case anisodamagerot:
04081 case creepbaz:
04082 case creepeffym:
04083 case contmat:
04084 case cebfipcontmat:
04085 case consolidation:
04086 case nonlocplastmat:
04087 case nonlocdamgmat:
04088 case nonlocalmicroM4:
04089 case glasgowdamage:
04090 case relaxationeuro:
04091 case simvisc:
04092 case lemaitr:
04093 case glasgowmechmat:
04094 case viscoplasticity:
04095 case lenjonespot:
04096 case creep_damage:
04097 case creepdpl:
04098 case damage_plasticity:
04099 break;
04100
04101
04102
04103
04104 case creeprs:
04105 crrs[im].give_reqnmq(anmq);
04106 break;
04107 case creepb3:
04108 crb3[im].give_reqnmq(anmq);
04109 break;
04110 case modcamclaycoupmat:
04111 cclayc[im].give_reqnmq(anmq);
04112 break;
04113 case hypoplastmat:
04114 hypopl[im].give_reqnmq(anmq);
04115 break;
04116 case therisodilat:
04117 tidilat[im].give_reqnmq(anmq);
04118 break;
04119 case effective_stress:
04120 effstr[im].give_reqnmq(anmq);
04121 break;
04122 case shrinkagemat:
04123 shmat[im].give_reqnmq(anmq);
04124 break;
04125 case time_switchmat:
04126
04127 break;
04128
04129 default:
04130 print_err("unknown material type is required", __FILE__, __LINE__, __func__);
04131 }
04132
04133 return;
04134 }
04135
04136
04137
04138
04139
04140
04141
04142
04143
04144
04145
04146
04147
04148
04149
04150 double mechmat::givenonmechq(nonmechquant qt, long ipp)
04151 {
04152 long id = nmqid[qt-1];
04153 if (id < 0)
04154 {
04155 print_err("Required quantity %s is not used in the problem solved",
04156 __FILE__, __LINE__, __func__, nonmechquantstr[qt-1].alias);
04157 abort();
04158 }
04159 return nonmechq[id][ipp];
04160 }
04161
04162
04163
04164
04165
04166
04167
04168
04169
04170
04171
04172
04173
04174
04175
04176
04177 void mechmat::storenonmechq(nonmechquant qt, long ipp, double val)
04178 {
04179 long id = nmqid[qt-1];
04180 if (id < 0)
04181 {
04182 print_err("Required quantity %s is not used in the problem solved",
04183 __FILE__, __LINE__, __func__, nonmechquantstr[qt-1].alias);
04184 abort();
04185 }
04186 nonmechq[id][ipp] = val;
04187 }
04188
04189
04190
04191
04192
04193
04194
04195
04196
04197
04198
04199
04200
04201
04202 long mechmat::givestatusnmq (nonmechquant qt)
04203 {
04204 if (nmqid[qt-1] < 0)
04205 return 0;
04206
04207 return 1;
04208 }
04209
04210
04211
04212
04213
04214
04215
04216
04217
04218
04219
04220
04221
04222
04223 long mechmat::givenonmechqid (nonmechquant qt)
04224 {
04225 return nmqid[qt-1];
04226 }
04227
04228
04229
04230
04231
04232
04233
04234
04235
04236
04237
04238
04239
04240
04241
04242
04243
04244
04245
04246
04247
04248
04249
04250
04251
04252
04253
04254
04255
04256 void mechmat::matstiff (matrix &d,long ipp,long im,long ido)
04257 {
04258 long i;
04259 switch (ip[ipp].tm[im]){
04260 case elisomat:{
04261 i=ip[ipp].idm[im];
04262 eliso[i].matstiff (d,ip[ipp].ssst);
04263 break;
04264 }
04265 case elgmat3d:{
04266 i=ip[ipp].idm[im];
04267 elgm3d[i].matstiff (d);
04268 break;
04269 }
04270 case simplas1d:{
04271 i=ip[ipp].idm[im];
04272 spl1d[i].matstiff (d,ipp,ido);
04273 break;
04274 }
04275 case jflow:{
04276 i=ip[ipp].idm[im];
04277 j2f[i].matstiff (d,ipp,ido);
04278 break;
04279 }
04280 case jflow2:{
04281 i=ip[ipp].idm[im];
04282 j2f2[i].matstiff (d,ipp,ido);
04283 break;
04284 }
04285 case microplaneM4:{
04286 i=ip[ipp].idm[im];
04287 mpM4[i].matstiff (d,ipp,ido);
04288 break;
04289 }
04290 case microsimp:{
04291 i=ip[ipp].idm[im];
04292 mpSIM[i].matstiff (d,ipp,ido);
04293 break;
04294 }
04295 case microfibro:{
04296 i=ip[ipp].idm[im];
04297 mpfib[i].matstiff (d,ipp,ido);
04298 break;
04299 }
04300 case mohrcoul:{
04301 i=ip[ipp].idm[im];
04302 mcoul[i].matstiff (d,ipp,ido);
04303 break;
04304 }
04305 case mohrcoulparab:{
04306 i=ip[ipp].idm[im];
04307 mcpar[i].matstiff (d,ipp,ido);
04308 break;
04309 }
04310 case boermaterial:{
04311 i=ip[ipp].idm[im];
04312 boerm[i].matstiff (d,ipp,ido);
04313 break;
04314 }
04315 case druckerprager:{
04316 i=ip[ipp].idm[im];
04317 drprm[i].matstiff (d,ipp,ido);
04318 break;
04319 }
04320 case doubledrprager:{
04321 i=ip[ipp].idm[im];
04322 ddpm[i].matstiff (d,ipp,ido);
04323 break;
04324 }
04325 case druckerprager2:{
04326 i=ip[ipp].idm[im];
04327 drprm2[i].matstiff (d,ipp,ido);
04328 break;
04329 }
04330 case chenplast:{
04331 i=ip[ipp].idm[im];
04332 chplast[i].matstiff (d,ipp,ido);
04333 break;
04334 }
04335 case modcamclaymat:{
04336 i=ip[ipp].idm[im];
04337 cclay[i].matstiff (d,ipp,im,ido);
04338 break;
04339 }
04340 case modcamclaycoupmat:{
04341 i=ip[ipp].idm[im];
04342 cclayc[i].matstiff (d,ipp,ido);
04343 break;
04344 }
04345 case hypoplastmat:{
04346 i=ip[ipp].idm[im];
04347 hypopl[i].matstiff (d,ipp,ido);
04348 break;
04349 }
04350 case shefpl:{
04351 i=ip[ipp].idm[im];
04352 shpl[i].matstiff (d,ipp,ido);
04353 break;
04354 }
04355 case hissplasmat:{
04356 i=ip[ipp].idm[im];
04357 hisspl[i].matstiff (d,ipp,ido);
04358 break;
04359 }
04360 case glasgowmechmat:{
04361 i=ip[ipp].idm[im];
04362 glasgmat[i].matstiff (d,ipp,ido);
04363 break;
04364 }
04365 case glasgowdamage:{
04366 i=ip[ipp].idm[im];
04367 glasgdam[i].matstiff (d,ipp,ido);
04368 break;
04369 }
04370 case creep_damage:{
04371 i=ip[ipp].idm[im];
04372 crdam[i].matstiff (d,ipp,im,ido);
04373 break;
04374 }
04375 case time_switchmat:{
04376 i=ip[ipp].idm[im];
04377 tswmat[i].matstiff (d,ipp,im,ido);
04378 break;
04379 }
04380 case effective_stress:{
04381 matstiff (d,ipp,im+1,ido);
04382 break;
04383 }
04384
04385
04386
04387
04388
04389
04390
04391 case scaldamage:{
04392 i=ip[ipp].idm[im];
04393 scdam[i].matstiff (d,ipp,ido);
04394 break;
04395 }
04396 case fixortodamage:{
04397 i=ip[ipp].idm[im];
04398 fixodam[i].matstiff (d,ipp,ido);
04399 break;
04400 }
04401 case scaldamagecc:{
04402 i=ip[ipp].idm[im];
04403 scdamcc[i].matstiff (d,ipp,ido);
04404 break;
04405 }
04406 case ortodamage:{
04407 i=ip[ipp].idm[im];
04408 ortdam[i].matstiff (d,ipp,ido);
04409 break;
04410 }
04411 case ortodamagerot:{
04412 i=ip[ipp].idm[im];
04413 ortdamrot[i].matstiff (d,ipp,ido);
04414 break;
04415 }
04416 case anisodamage:{
04417 i=ip[ipp].idm[im];
04418 anidam[i].matstiff (d,ipp,ido);
04419 break;
04420 }
04421 case anisodamagerot:{
04422 i=ip[ipp].idm[im];
04423 anidamrot[i].matstiff (d,ipp,ido);
04424 break;
04425 }
04426 case graphm:{
04427 i=ip[ipp].idm[im];
04428 grmat[i].matstiff (d,ipp);
04429 break;
04430 }
04431 case varelisomat:{
04432 i=ip[ipp].idm[im];
04433 veliso[i].matstiff (d,ipp);
04434 break;
04435 }
04436 case geoelast:{
04437 i=ip[ipp].idm[im];
04438 geoel[i].matstiff (d,ipp,ido);
04439 break;
04440 }
04441 case creeprs:
04442 case creepb3:{
04443 creep_matstiff (d,ipp,im,ido);
04444 break;
04445 }
04446 case creepdpl:{
04447 creep_matstiff (d,ipp,im,ido);
04448 break;
04449 }
04450 case creepbaz:{
04451 i=ip[ipp].idm[im];
04452 crbaz[i].matstiff (d,ipp);
04453 break;
04454 }
04455 case creepeffym:{
04456 i=ip[ipp].idm[im];
04457 creffym[i].matstiff (d,ipp,im,ido);
04458 break;
04459 }
04460 case consolidation:{
04461 i=ip[ipp].idm[im];
04462 csol[i].matstiff (d,ipp);
04463 break;
04464 }
04465 case winklerpasternak:{
04466 i=ip[ipp].idm[im];
04467 wpast[i].matstiff (d,ipp);
04468 break;
04469 }
04470
04471 case nonlocdamgmat:{
04472 matstiff (d,ipp,im+1,ido);
04473 break;
04474 }
04475 case nonlocplastmat:{
04476 matstiff (d,ipp,im+1,ido);
04477 break;
04478 }
04479 case nonlocalmicroM4:{
04480 matstiff (d,ipp,im+1,ido);
04481 break;
04482 }
04483
04484 case damage_plasticity:{
04485 i=ip[ipp].idm[im];
04486 dampl[i].matstiff (d,ipp);
04487 break;
04488 }
04489
04490 case viscoplasticity:{
04491 i=ip[ipp].idm[im];
04492 visplas[i].matstiff (d,ipp,im,ido);
04493 break;
04494 }
04495
04496 case contmat:{
04497 i=ip[ipp].idm[im];
04498 conmat[i].matstiff (d,ipp);
04499 break;
04500 }
04501
04502 case cebfipcontmat:{
04503 i=ip[ipp].idm[im];
04504 cebfipconmat[i].matstiff (d,ipp);
04505 break;
04506 }
04507 case shrinkagemat:{
04508 i=ip[ipp].idm[im];
04509 shmat[i].matstiff (d,ipp,im,ido);
04510 break;
04511 }
04512
04513 default:
04514 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
04515
04516 }
04517 }
04518
04519
04520
04521
04522
04523
04524
04525
04526
04527
04528
04529
04530
04531
04532
04533 void mechmat::elmatstiff (matrix &d,long ipp)
04534
04535 {
04536 long i,idem;
04537
04538 idem = ip[ipp].gemid();
04539 switch (ip[ipp].tm[idem]){
04540 case elisomat:{
04541 i=ip[ipp].idm[idem];
04542 eliso[i].elmatstiff (d,ip[ipp].ssst);
04543 break;
04544 }
04545 case elortomat:{
04546 i=ip[ipp].idm[idem];
04547 elorto[i].elmatstiff (d,ip[ipp].ssst);
04548 break;
04549 }
04550 case elgmat3d:{
04551 i=ip[ipp].idm[idem];
04552 elgm3d[i].elmatstiff (d);
04553 break;
04554 }
04555 case varelisomat:{
04556 i=ip[ipp].idm[idem];
04557 veliso[i].elmatstiff (d,ipp);
04558 break;
04559 }
04560 default:
04561 print_err("Unknown elastic material type is required",__FILE__,__LINE__,__func__);
04562 }
04563 }
04564
04565
04566
04567
04568
04569
04570
04571
04572
04573
04574
04575
04576
04577
04578 void mechmat::elmatcompl (matrix &c,long ipp)
04579
04580 {
04581 long i, idem;
04582
04583 idem = ip[ipp].gemid();
04584 switch (ip[ipp].tm[idem]){
04585 case elisomat:{
04586 i=ip[ipp].idm[idem];
04587 eliso[i].matcompl (c,ip[ipp].ssst);
04588 break;
04589 }
04590 case varelisomat:{
04591 i=ip[ipp].idm[idem];
04592 veliso[i].matcompl (c,ipp);
04593 break;
04594 }
04595 default:
04596 print_err("unknown elastic material type is required",__FILE__,__LINE__,__func__);
04597
04598 }
04599 }
04600
04601
04602
04603
04604
04605
04606
04607
04608
04609
04610
04611
04612
04613
04614 void mechmat::eigenstresses (vector &sig,long ipp,long id)
04615 {
04616 long i,n;
04617 n=sig.n;
04618 vector eps;
04619 allocv (n,eps);
04620
04621 switch (ip[ipp].tm[id]){
04622 case relaxationeuro:{
04623
04624 i=ip[ipp].idm[id];
04625
04626 relaxec[i].stress (sig,eps,ip[ipp].ssst);
04627
04628 storeeigstrain (ipp,eps);
04629
04630 storeeigstress (ipp,sig);
04631 break;
04632 }
04633 default:
04634 print_err("Unknown model for eigenstresses computation is required",__FILE__,__LINE__,__func__);
04635 }
04636 destrv (eps);
04637 }
04638
04639
04640
04641
04642
04643
04644
04645
04646
04647
04648
04649
04650
04651
04652
04653
04654 void mechmat::computenlstresses (long ipp,long im,long ido)
04655 {
04656 mattype tm = ip[ipp].tm[im];
04657 long idm = ip[ipp].idm[im];
04658 computenlstressesmat(tm, idm, ipp, im, ido);
04659 }
04660
04661
04662 void mechmat::computenlstressesmat (mattype tm, long idm, long ipp,long im,long ido)
04663 {
04664 switch (tm){
04665 case elisomat:{
04666 eliso[idm].nlstresses (ipp);
04667 break;
04668 }
04669 case elgmat3d:{
04670 elgm3d[ip[ipp].idm[im]].nlstresses (ipp);
04671 break;
04672 }
04673 case simplas1d:{
04674 spl1d[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04675 break;
04676 }
04677 case jflow:{
04678 j2f[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04679 break;
04680 }
04681 case jflow2:{
04682 j2f2[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04683 break;
04684 }
04685 case microplaneM4:{
04686 mpM4[ip[ipp].idm[im]].nlstresses (ipp,ido);
04687 break;
04688 }
04689 case microsimp:{
04690 mpSIM[ip[ipp].idm[im]].nlstresses (ipp,ido);
04691 break;
04692 }
04693 case microfibro:{
04694 mpfib[ip[ipp].idm[im]].nlstresses (ipp,ido);
04695 break;
04696 }
04697 case mohrcoul:{
04698 mcoul[ip[ipp].idm[im]].nlstresses (ipp,ido);
04699 break;
04700 }
04701 case mohrcoulparab:{
04702 mcpar[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04703 break;
04704 }
04705 case boermaterial:{
04706 boerm[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04707 break;
04708 }
04709 case druckerprager:{
04710 drprm[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04711 break;
04712 }
04713 case doubledrprager:{
04714 ddpm[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04715 break;
04716 }
04717 case druckerprager2:{
04718 drprm2[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04719 break;
04720 }
04721 case chenplast:{
04722 chplast[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04723 break;
04724 }
04725 case modcamclaymat:{
04726 cclay[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04727 break;
04728 }
04729 case modcamclaycoupmat:{
04730 cclayc[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04731 break;
04732 }
04733 case hypoplastmat:{
04734 hypopl[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04735 break;
04736 }
04737 case hissplasmat:{
04738 hisspl[ip[ipp].idm[im]].nlstresses (ipp,ido);
04739 break;
04740 }
04741 case shefpl:{
04742 shpl[ip[ipp].idm[im]].nlstresses (ipp,ido);
04743 break;
04744 }
04745 case glasgowmechmat:{
04746 glasgmat[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04747 break;
04748 }
04749 case glasgowdamage:{
04750 glasgdam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04751 break;
04752 }
04753 case creep_damage:{
04754 crdam[ip[ipp].idm[im]].nlstresses (ipp, im, ido);
04755 break;
04756 }
04757 case shrinkagemat:{
04758 shmat[ip[ipp].idm[im]].nlstresses (ipp, im, ido);
04759 break;
04760 }
04761 case time_switchmat:{
04762 tswmat[ip[ipp].idm[im]].nlstresses (ipp, im, ido);
04763 break;
04764 }
04765 case effective_stress:{
04766 effstr[ip[ipp].idm[im]].nlstresses (ipp, im, ido);
04767 break;
04768 }
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779
04780
04781 case scaldamage:{
04782 scdam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04783 break;
04784 }
04785 case fixortodamage:{
04786 fixodam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04787 break;
04788 }
04789 case scaldamagecc:{
04790 scdamcc[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04791 break;
04792 }
04793 case ortodamage:{
04794 ortdam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04795 break;
04796 }
04797 case ortodamagerot:{
04798 ortdamrot[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04799 break;
04800 }
04801 case anisodamage:{
04802 anidam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04803 break;
04804 }
04805 case anisodamagerot:{
04806 anidamrot[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04807 break;
04808 }
04809 case graphm:{
04810 grmat[ip[ipp].idm[im]].nlstresses (ipp);
04811 break;
04812 }
04813 case varelisomat:{
04814 veliso[ip[ipp].idm[im]].nlstresses (ipp);
04815 break;
04816 }
04817 case geoelast:{
04818 geoel[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04819 break;
04820 }
04821 case creepb3:
04822 case creeprs:
04823 case creepdpl:{
04824 creep_nlstresses(ipp,im,ido);
04825 break;
04826 }
04827 case creepbaz:{
04828 crbaz[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04829 break;
04830 }
04831 case creepeffym:{
04832 creffym[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04833 break;
04834 }
04835 case consolidation:{
04836 csol[ip[ipp].idm[im]].nlstresses (ipp);
04837 break;
04838 }
04839 case winklerpasternak:{
04840 wpast[ip[ipp].idm[im]].nlstresses (ipp);
04841 break;
04842 }
04843
04844 case nonlocdamgmat:{
04845 computenlstresses (ipp,im+1,ido);
04846 break;
04847 }
04848 case nonlocplastmat:{
04849 computenlstresses (ipp,im+1,ido);
04850 break;
04851 }
04852 case nonlocalmicroM4:{
04853 computenlstresses (ipp,im+1,ido);
04854 break;
04855 }
04856 case damage_plasticity:{
04857 dampl[ip[ipp].idm[im]].nlstresses (ipp);
04858 break;
04859 }
04860 case viscoplasticity:{
04861 visplas[ip[ipp].idm[im]].nlstresses (ipp,ido);
04862 break;
04863 }
04864
04865
04866 case contmat:{
04867 conmat[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04868 break;
04869 }
04870
04871 case cebfipcontmat:{
04872 cebfipconmat[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
04873 break;
04874 }
04875
04876 default:{
04877 print_err("unknown material type is required", __FILE__,__LINE__,__func__);
04878 }
04879 }
04880
04881 }
04882
04883
04884
04885
04886
04887
04888
04889
04890
04891
04892
04893
04894
04895
04896
04897 void mechmat::computenlstressesincr (long ipp,long im,long ido)
04898 {
04899 switch (ip[ipp].tm[im]){
04900
04901 case creep_damage:{
04902 crdam[ip[ipp].idm[im]].nlstressesincr (ipp, im, ido);
04903 break;
04904 }
04905 case shrinkagemat:{
04906 shmat[ip[ipp].idm[im]].nlstressesincr (ipp, im, ido);
04907 break;
04908 }
04909
04910
04911
04912
04913
04914
04915
04916
04917
04918
04919
04920
04921 case creepb3:
04922 case creeprs:
04923 case creepdpl:{
04924 creep_nlstressesincr (ipp,im,ido);
04925 break;
04926 }
04927
04928
04929
04930
04931
04932
04933
04934
04935
04936
04937
04938 case consolidation:{
04939 csol[ip[ipp].idm[im]].nlstressesincr (ipp);
04940 break;
04941 }
04942
04943
04944 case viscoplasticity:{
04945 visplas[ip[ipp].idm[im]].nlstressesincr (ipp,im,ido);
04946 break;
04947 }
04948
04949 case elisomat:
04950 case simplas1d:
04951 case jflow:
04952 case jflow2:
04953 case microplaneM4:
04954 case microsimp:
04955 case microfibro:
04956 case mohrcoul:
04957 case mohrcoulparab:
04958 case boermaterial:
04959 case druckerprager:
04960 case doubledrprager:
04961 case druckerprager2:
04962 case chenplast:
04963 case modcamclaymat:
04964 case modcamclaycoupmat:
04965 case hypoplastmat:
04966 case shefpl:
04967 case effective_stress:
04968 case glasgowmechmat:
04969 case glasgowdamage:
04970 case scaldamage:
04971 case fixortodamage:
04972 case scaldamagecc:
04973 case ortodamage:
04974 case ortodamagerot:
04975 case anisodamage:
04976 case anisodamagerot:
04977 case graphm:
04978 case varelisomat:
04979 case geoelast:
04980 case creepeffym:
04981 case winklerpasternak:
04982 case nonlocdamgmat:
04983 case nonlocplastmat:
04984 case nonlocalmicroM4:
04985 case damage_plasticity:
04986 case cebfipcontmat:
04987 case contmat:
04988 if (im == 0)
04989 ip[ipp].clean_stresses(Mb->nlc);
04990 break;
04991 default:{
04992 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
04993 }
04994 }
04995
04996 }
04997
04998
04999
05000
05001
05002
05003
05004
05005
05006
05007
05008
05009
05010
05011
05012 void mechmat::compnonloc_nlstresses (long ipp,long im,long ido)
05013 {
05014 switch (ip[ipp].tm[im]){
05015 case elisomat:{
05016 eliso[ip[ipp].idm[im]].nlstresses (ipp);
05017 break;
05018 }
05019 case simplas1d:{
05020 if (ip[ipp].hmt & 2)
05021 spl1d[ip[ipp].idm[im]].nonloc_nlstresses (ipp,im,ido);
05022 else
05023 spl1d[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05024 break;
05025 }
05026 case jflow:{
05027 if (ip[ipp].hmt & 2)
05028 j2f[ip[ipp].idm[im]].nonloc_nlstresses (ipp,im,ido);
05029 else
05030 j2f[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05031 break;
05032 }
05033 case microplaneM4:{
05034 if (ip[ipp].hmt & 2)
05035 mpM4[ip[ipp].idm[im]].nonloc_nlstresses (ipp,ido);
05036 else
05037 mpM4[ip[ipp].idm[im]].nlstresses (ipp,ido);
05038 break;
05039 }
05040 case microsimp:{
05041
05042 break;
05043 }
05044 case microfibro:{
05045
05046 break;
05047 }
05048 case mohrcoul:{
05049 if (ip[ipp].hmt & 2)
05050 mcoul[ip[ipp].idm[im]].nonloc_nlstresses (ipp,ido);
05051 else
05052 mcoul[ip[ipp].idm[im]].nlstresses (ipp,ido);
05053 break;
05054 }
05055 case mohrcoulparab:{
05056 if (ip[ipp].hmt & 2)
05057 mcpar[ip[ipp].idm[im]].nonloc_nlstresses (ipp,im,ido);
05058 else
05059 mcpar[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05060 break;
05061 }
05062 case boermaterial:{
05063 if (ip[ipp].hmt & 2)
05064 boerm[ip[ipp].idm[im]].nonloc_nlstresses (ipp,im,ido);
05065 else
05066 boerm[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05067 break;
05068 }
05069 case druckerprager:{
05070 if (ip[ipp].hmt & 2)
05071 drprm[ip[ipp].idm[im]].nonloc_nlstresses (ipp,im,ido);
05072 else
05073 drprm[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05074 break;
05075 }
05076
05077 case doubledrprager:{
05078
05079
05080
05081
05082 break;
05083 }
05084 case druckerprager2:{
05085
05086
05087
05088
05089 break;
05090 }
05091 case modcamclaymat:{
05092
05093 break;
05094 }
05095 case hissplasmat:{
05096
05097 break;
05098 }
05099 case effective_stress:{
05100 compnonloc_nlstresses(ipp, im+1, ido);
05101 break;
05102 }
05103 case simvisplas:{
05104
05105 break;
05106 }
05107 case lemaitr:{
05108
05109 break;
05110 }
05111 case scaldamage:{
05112 scdam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05113 break;
05114 }
05115 case fixortodamage:{
05116 fixodam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05117 break;
05118 }
05119 case scaldamagecc:{
05120 scdamcc[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05121 break;
05122 }
05123 case ortodamage:{
05124 ortdam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05125 break;
05126 }
05127 case ortodamagerot:{
05128 ortdamrot[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05129 break;
05130 }
05131 case anisodamage:{
05132 anidam[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05133 break;
05134 }
05135 case anisodamagerot:{
05136 anidamrot[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05137 break;
05138 }
05139 case graphm:{
05140 grmat[ip[ipp].idm[im]].nlstresses (ipp);
05141 break;
05142 }
05143 case varelisomat:{
05144 veliso[ip[ipp].idm[im]].nlstresses (ipp);
05145 break;
05146 }
05147 case creepbaz:{
05148
05149 break;
05150 }
05151 case consolidation:{
05152 csol[ip[ipp].idm[im]].nlstresses (ipp);
05153 break;
05154 }
05155 case winklerpasternak:{
05156 wpast[ip[ipp].idm[im]].nlstresses (ipp);
05157 break;
05158 }
05159
05160 case creep_damage:{
05161 if (ip[ipp].hmt & 2)
05162 crdam[ip[ipp].idm[im]].nonloc_nlstresses (ipp, im, ido);
05163 else
05164 crdam[ip[ipp].idm[im]].nlstresses (ipp, im, ido);
05165 break;
05166 }
05167 case shrinkagemat:{
05168 if (ip[ipp].hmt & 2)
05169 shmat[ip[ipp].idm[im]].nonloc_nlstresses (ipp, im, ido);
05170 else
05171 shmat[ip[ipp].idm[im]].nlstresses (ipp, im, ido);
05172 break;
05173 }
05174 case creepb3:
05175 case creeprs:
05176 case creepdpl:{
05177 creep_nlstresses(ipp,im,ido);
05178 break;
05179 }
05180 case creepeffym:{
05181 creffym[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
05182 break;
05183 }
05184 case nonlocdamgmat:{
05185 compnonloc_nlstresses (ipp,im+1,ido);
05186 break;
05187 }
05188 case nonlocplastmat:{
05189 compnonloc_nlstresses (ipp,im+1,ido);
05190 break;
05191 }
05192 case nonlocalmicroM4:{
05193 compnonloc_nlstresses (ipp,im+1,ido);
05194 break;
05195 }
05196
05197 default:{
05198 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
05199 }
05200 }
05201
05202 }
05203
05204
05205
05206
05207
05208
05209
05210
05211
05212
05213
05214
05215 void mechmat::compute_averstrains ()
05216 {
05217 long i,ne;
05218 double volume;
05219 elemtype te;
05220 vector averstra (6);
05221
05222 volume=0.0;
05223 fillv (0.0,averstra);
05224
05225 ne=Mt->ne;
05226 for (i=0;i<ne;i++){
05227 if (Gtm->leso[i]==1){
05228
05229 te = Mt->give_elem_type (i);
05230
05231 switch (te){
05232 case linearhex:{
05233 Lhex->aver_strains (0,i,0,0,averstra,volume);
05234 break;
05235 }
05236 case quadrhex:{
05237 Qhex->aver_strains (0,i,0,0,averstra,volume);
05238 break;
05239 }
05240
05241 default:{
05242
05243 }
05244 }
05245 }
05246 }
05247
05248 for (i=0;i<6;i++){
05249 averstra[i]/=volume;
05250 }
05251
05252 printf ("\n\n prumerne deformace");
05253 for (i=0;i<6;i++){
05254 printf ("\n slozka %ld %f",i+1,averstra[i]);
05255 }
05256
05257 }
05258
05259
05260
05261
05262
05263
05264
05265
05266
05267
05268
05269
05270 void mechmat::temprstrains (long lcid)
05271 {
05272 long i,nm,ii;
05273
05274 for (i=0;i<tnip;i++){
05275
05276 ii = ip[i].hmt & 1;
05277 if (ii > 0){
05278 nm=ip[i].nm-1;
05279 switch (ip[i].tm[nm]){
05280 case therisodilat:{
05281 tidilat[ip[i].idm[nm]].temprstrains (lcid,i);
05282 break;
05283 }
05284 default:
05285 print_err("unknown thermal dilatation material is required",__FILE__,__LINE__,__func__);
05286 }
05287 }
05288
05289 }
05290 }
05291
05292
05293
05294
05295
05296
05297
05298
05299
05300
05301 void mechmat::eigstrmod ()
05302 {
05303 long i,j,k,nc,ipp,nip;
05304
05305 if (eigstrains != NULL){
05306 for (i=0;i<Mt->ne;i++){
05307 if (Gtm->leso[i]==1){
05308
05309 ipp=Mt->elements[i].ipp[0][0];
05310
05311 nip = Mt->give_tnip (i);
05312 for (j=0;j<nip;j++){
05313
05314 nc=ip[ipp].ncompstr;
05315 for (k=0;k<nc;k++)
05316 ip[ipp].strain[k]-=eigstrains[ipp][k];
05317 ipp++;
05318 }
05319 }
05320 }
05321 }
05322 if (tempstrains != NULL){
05323 for (i=0;i<Mt->ne;i++){
05324 if (Gtm->leso[i]==1){
05325
05326 ipp=Mt->elements[i].ipp[0][0];
05327
05328 nip = Mt->give_tnip (i);
05329 for (j=0;j<nip;j++){
05330
05331 nc=ip[ipp].ncompstr;
05332 for (k=0;k<nc;k++)
05333 ip[ipp].strain[k]-=tempstrains[ipp][k];
05334 ipp++;
05335 }
05336 }
05337 }
05338 }
05339 }
05340
05341
05342
05343
05344
05345
05346
05347
05348
05349
05350
05351 void mechmat::add_macro_strains (vector ¯ostrains)
05352 {
05353 long i,j,ncompstr;
05354
05355 for (i=0;i<tnip;i++){
05356 ncompstr = ip[i].ncompstr;
05357 if (ncompstr != macrostrains.n){
05358 print_err("the number of strain components (%ld) in integration point %ld\n"
05359 "is not equal to the number of macro-strain component (%ld) in homogenization",
05360 __FILE__, __LINE__, __func__,ncompstr,i,macrostrains.n);
05361 abort ();
05362 }
05363 for (j=0;j<ncompstr;j++){
05364 ip[i].strain[j] += macrostrains[j];
05365 }
05366 }
05367 }
05368
05369
05370
05371
05372
05373
05374
05375
05376 void mechmat::totstrains ()
05377 {
05378 long i,j,k,nc,ipp,nip;
05379
05380 if (eigstrains != NULL){
05381 for (i=0;i<Mt->ne;i++){
05382 if (Gtm->leso[i]==1){
05383
05384 ipp=Mt->elements[i].ipp[0][0];
05385
05386 nip = Mt->give_tnip (i);
05387 for (j=0;j<nip;j++){
05388
05389 nc=ip[ipp].ncompstr;
05390 for (k=0;k<nc;k++)
05391 ip[ipp].strain[k]+=eigstrains[ipp][k];
05392 ipp++;
05393 }
05394 }
05395 }
05396 }
05397 if (tempstrains != NULL){
05398 for (i=0;i<Mt->ne;i++){
05399 if (Gtm->leso[i]==1){
05400
05401 ipp=Mt->elements[i].ipp[0][0];
05402
05403 nip = Mt->give_tnip (i);
05404 for (j=0;j<nip;j++){
05405
05406 nc=ip[ipp].ncompstr;
05407 for (k=0;k<nc;k++)
05408 ip[ipp].strain[k]+=tempstrains[ipp][k];
05409 ipp++;
05410 }
05411 }
05412 }
05413 }
05414 }
05415
05416
05417
05418
05419
05420
05421
05422
05423
05424
05425
05426
05427
05428
05429 void mechmat::matdilat (matrix &d,long ipp)
05430 {
05431 long i, iddm = ip[ipp].nm-1;
05432
05433 switch (ip[ipp].tm[iddm]){
05434 case therisodilat:{
05435 i=ip[ipp].idm[iddm];
05436 tidilat[i].matdilat (d,ip[ipp].ssst);
05437 break;
05438 }
05439 default:{
05440 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
05441 }
05442 }
05443
05444 }
05445
05446
05447
05448
05449
05450
05451
05452
05453
05454
05455
05456
05457
05458 void mechmat::storeipvol (long ipp,double vol)
05459 {
05460 ipv[ipp]=vol;
05461 }
05462
05463
05464
05465
05466
05467
05468
05469
05470
05471
05472
05473
05474
05475
05476 void mechmat::nonlocaverage (long ipp,long im,long ido)
05477 {
05478 long tncompo;
05479 long nldamncompo;
05480
05481 if (ip[ipp].hmt & 2)
05482 {
05483 switch (ip[ipp].tm[im]){
05484 case nonlocplastmat :
05485 nlplast[ip[ipp].idm[im]].average(ipp, ido);
05486 break;
05487 case nonlocdamgmat :
05488 nldamg[ip[ipp].idm[im]].average(ipp,im,ido);
05489 break;
05490 case nonlocalmicroM4 :
05491 nmpM4[ip[ipp].idm[im]].average(ipp);
05492 break;
05493 case creep_damage:
05494 tncompo = givencompeqother(ipp,im);
05495 nldamncompo = givencompeqother(ipp,im+2);
05496 nonlocaverage(ipp, im+2, ido+tncompo-nldamncompo);
05497 break;
05498 case shrinkagemat:
05499 tncompo = givencompeqother(ipp,im);
05500 nldamncompo = givencompeqother(ipp,im+1);
05501 nonlocaverage(ipp, im+1, ido+tncompo-nldamncompo);
05502 break;
05503 default:{
05504 print_err("unknown type of averaging model is required", __FILE__, __LINE__, __func__);
05505 }
05506 }
05507 }
05508 }
05509
05510
05511
05512
05513
05514
05515
05516
05517
05518
05519
05520
05521
05522
05523 long mechmat::give_num_averq(long ipp,long im)
05524 {
05525 long ret = 0;
05526 if (ip[ipp].hmt & 2)
05527 {
05528 switch (ip[ipp].tm[im]){
05529 case nonlocalmicroM4 :
05530 ret = nmpM4[ip[ipp].idm[im]].give_num_averq(ipp);
05531 break;
05532 case nonlocplastmat :
05533 ret = nlplast[ip[ipp].idm[im]].give_num_averq(ipp);
05534 break;
05535 case nonlocdamgmat :
05536 ret = nldamg[ip[ipp].idm[im]].give_num_averq(ipp, im);
05537 break;
05538 case creep_damage:
05539 ret = give_num_averq(ipp, im+2);
05540 break;
05541 case shrinkagemat:
05542 ret = give_num_averq(ipp, im+1);
05543 break;
05544 default:{
05545 print_err("unknown type of averaging model is required", __FILE__, __LINE__, __func__);
05546 }
05547 }
05548 }
05549 return ret;
05550 }
05551
05552
05553
05554
05555
05556
05557
05558
05559
05560
05561
05562
05563
05564
05565 void mechmat::give_aver_quantv(long ipp,long im,long ido, vector &qv)
05566 {
05567 long tncompo;
05568 long nldamncompo;
05569
05570 if (ip[ipp].hmt & 2)
05571 {
05572 switch (ip[ipp].tm[im]){
05573 case nonlocplastmat :
05574 nlplast[ip[ipp].idm[im]].give_aver_quantv(ipp, im, ido, qv);
05575 break;
05576 case nonlocdamgmat :
05577 nldamg[ip[ipp].idm[im]].give_aver_quantv(ipp, im, ido, qv);
05578 break;
05579 case creep_damage:
05580 tncompo = givencompeqother(ipp,im);
05581 nldamncompo = givencompeqother(ipp,im+2);
05582 give_aver_quantv(ipp, im+2, ido+tncompo-nldamncompo, qv);
05583 break;
05584 case shrinkagemat:
05585 tncompo = givencompeqother(ipp,im);
05586 nldamncompo = givencompeqother(ipp,im+2);
05587 give_aver_quantv(ipp, im+2, ido+tncompo-nldamncompo, qv);
05588 break;
05589 default:{
05590 print_err("unknown type of averaging model is required", __FILE__, __LINE__, __func__);
05591 }
05592 }
05593 }
05594 return;
05595 }
05596
05597
05598
05599
05600
05601
05602
05603
05604
05605
05606
05607
05608
05609 double mechmat::nonlocradius (long ipp,long im)
05610 {
05611 double r;
05612
05613 switch (ip[ipp].tm[im]){
05614 case nonlocplastmat :
05615 {
05616 r = nlplast[ip[ipp].idm[im]].r;
05617 break;
05618 }
05619 case nonlocdamgmat :
05620 {
05621 r = nldamg[ip[ipp].idm[im]].r;
05622 break;
05623 }
05624 case nonlocalmicroM4 :
05625 {
05626 r = nmpM4[ip[ipp].idm[im]].r;
05627 break;
05628 }
05629 default:{
05630 print_err("wrong material type is required",__FILE__,__LINE__,__func__);
05631 }
05632 }
05633
05634 return r;
05635 }
05636
05637
05638
05639
05640
05641
05642
05643
05644
05645
05646
05647
05648
05649 long mechmat::givenonlocid(long ipp)
05650 {
05651 long i, ret=-1;
05652 for(i=0; i<ip[ipp].nm; i++)
05653 {
05654 switch(ip[ipp].tm[i])
05655 {
05656 case nonlocplastmat:
05657 case nonlocdamgmat:
05658 case nonlocalmicroM4:
05659 ret = i;
05660 break;
05661 default:
05662 break;
05663 }
05664 if (ret >= 0)
05665 break;
05666 }
05667 return ret;
05668 }
05669
05670
05671
05672
05673
05674
05675
05676
05677
05678
05679
05680
05681 void mechmat::updateipval (void)
05682 {
05683 long i,j,ipp,nip;
05684
05685
05686 for (i=0;i<Mt->ne;i++){
05687 if (Gtm->leso[i]==1){
05688 ipp=Mt->elements[i].ipp[0][0];
05689 nip=Mt->give_tnip (i);
05690 for (j=0;j<nip;j++){
05691 updateipvalmat (ipp,0,0);
05692 ipp++;
05693 }
05694 }
05695 }
05696 }
05697
05698
05699
05700
05701
05702
05703
05704
05705
05706
05707
05708
05709
05710
05711 void mechmat::updateipvalmat (long ipp,long im,long ido)
05712 {
05713 switch (ip[ipp].tm[im]){
05714 case elisomat :
05715 break;
05716 case elortomat :
05717 break;
05718 case simplas1d:{
05719 spl1d[ip[ipp].idm[im]].updateval (ipp,ido);
05720 break;
05721 }
05722 case jflow:{
05723 j2f[ip[ipp].idm[im]].updateval (ipp,im,ido);
05724 break;
05725 }
05726 case jflow2:{
05727 j2f2[ip[ipp].idm[im]].updateval (ipp,im,ido);
05728 break;
05729 }
05730 case mohrcoul:{
05731 mcoul[ip[ipp].idm[im]].updateval (ipp,im,ido);
05732 break;
05733 }
05734 case mohrcoulparab:{
05735 mcpar[ip[ipp].idm[im]].updateval (ipp,im,ido);
05736 break;
05737 }
05738 case boermaterial:{
05739 boerm[ip[ipp].idm[im]].updateval (ipp,im,ido);
05740 break;
05741 }
05742 case druckerprager:{
05743 drprm[ip[ipp].idm[im]].updateval (ipp,im,ido);
05744 break;
05745 }
05746 case doubledrprager:{
05747 ddpm[ip[ipp].idm[im]].updateval (ipp,im,ido);
05748 break;
05749 }
05750 case druckerprager2:{
05751 drprm2[ip[ipp].idm[im]].updateval (ipp,im,ido);
05752 break;
05753 }
05754 case chenplast:{
05755 chplast[ip[ipp].idm[im]].updateval (ipp,im,ido);
05756 break;
05757 }
05758 case modcamclaymat:{
05759 cclay[ip[ipp].idm[im]].updateval (ipp,ido);
05760 break;
05761 }
05762 case modcamclaycoupmat:{
05763 cclayc[ip[ipp].idm[im]].updateval (ipp,ido);
05764 break;
05765 }
05766 case hypoplastmat:{
05767 hypopl[ip[ipp].idm[im]].updateval (ipp,im,ido);
05768 break;
05769 }
05770 case hissplasmat:{
05771 hisspl[ip[ipp].idm[im]].updateval (ipp,ido);
05772 break;
05773 }
05774 case effective_stress:
05775 updateipvalmat(ipp, im+1, ido);
05776 break;
05777 case microplaneM4:{
05778 mpM4[ip[ipp].idm[im]].updateval (ipp,im,ido);
05779 }
05780 case microsimp:{
05781 mpSIM[ip[ipp].idm[im]].updateval (ipp,im,ido);
05782 }
05783 case microfibro:{
05784 mpfib[ip[ipp].idm[im]].updateval (ipp,im,ido);
05785 }
05786 case shefpl:{
05787 shpl[ip[ipp].idm[im]].updateval (ipp,ido);
05788 break;
05789 }
05790 case glasgowdamage:{
05791 glasgdam[ip[ipp].idm[im]].updateval (ipp,im,ido);
05792 break;
05793 }
05794 case creep_damage:{
05795 crdam[ip[ipp].idm[im]].updateval (ipp, im, ido);
05796 break;
05797 }
05798 case shrinkagemat:{
05799 shmat[ip[ipp].idm[im]].updateval (ipp, im, ido);
05800 break;
05801 }
05802 case time_switchmat:{
05803 tswmat[ip[ipp].idm[im]].updateval (ipp, im, ido);
05804 break;
05805 }
05806 case graphm:{
05807 break;
05808 }
05809 case varelisomat :
05810 break;
05811 case geoelast:{
05812 geoel[ip[ipp].idm[im]].updateval (ipp,im,ido);
05813 break;
05814 }
05815 case scaldamage:{
05816 scdam[ip[ipp].idm[im]].updateval (ipp,im,ido);
05817 break;
05818 }
05819 case fixortodamage:{
05820 fixodam[ip[ipp].idm[im]].updateval (ipp,im,ido);
05821 break;
05822 }
05823 case scaldamagecc:{
05824 scdamcc[ip[ipp].idm[im]].updateval (ipp,im,ido);
05825 break;
05826 }
05827 case ortodamage:{
05828 ortdam[ip[ipp].idm[im]].updateval (ipp,im,ido);
05829 break;
05830 }
05831 case ortodamagerot:{
05832 ortdamrot[ip[ipp].idm[im]].updateval (ipp,im,ido);
05833 break;
05834 }
05835 case anisodamage:{
05836 anidam[ip[ipp].idm[im]].updateval (ipp,im,ido);
05837 break;
05838 }
05839 case anisodamagerot:{
05840 anidamrot[ip[ipp].idm[im]].updateval (ipp,im,ido);
05841 break;
05842 }
05843
05844 case nonlocdamgmat:{
05845 updateipvalmat (ipp,im+1,ido);
05846 break;
05847 }
05848 case nonlocplastmat:{
05849 updateipvalmat (ipp,im+1,ido);
05850 break;
05851 }
05852 case nonlocalmicroM4:{
05853 updateipvalmat (ipp,im+1,ido);
05854 break;
05855 }
05856
05857 case damage_plasticity:{
05858 dampl[ip[ipp].idm[im]].updateval (ipp);
05859 break;
05860 }
05861 case simvisc:
05862
05863 case lemaitr:
05864 case viscoplasticity:
05865 visplas[ip[ipp].idm[im]].updateval (ipp,im,ido);
05866 break;
05867 case creeprs:
05868 case creepb3:
05869 case creepdpl:{
05870 creep_updateval (ipp,im,ido);
05871 break;
05872 }
05873 case creepbaz:
05874 crbaz[ip[ipp].idm[im]].updateval ();
05875 break;
05876 case creepeffym:
05877 creffym[ip[ipp].idm[im]].updateval (ipp,im,ido);
05878 break;
05879 case lenjonespot:
05880 break;
05881
05882 case contmat:
05883 conmat[ip[ipp].idm[im]].updateval (ipp,im,ido);
05884 break;
05885
05886 case cebfipcontmat:
05887 cebfipconmat[ip[ipp].idm[im]].updateval (ipp,im,ido);
05888 break;
05889
05890
05891 case consolidation:{
05892 csol[ip[ipp].idm[im]].updateval (ipp,im,ido);
05893 break;
05894 }
05895 default:{
05896 print_err ("unknown material type is required",__FILE__,__LINE__,__func__);
05897 break;
05898 }
05899 }
05900 }
05901
05902
05903
05904
05905
05906
05907
05908
05909
05910
05911 void mechmat::clean_ip ()
05912 {
05913 long i;
05914 for (i=0;i<tnip;i++){
05915 ip[i].clean (Mb->nlc);
05916 }
05917 }
05918
05919
05920
05921
05922
05923
05924
05925
05926
05927
05928
05929
05930
05931
05932
05933
05934 double mechmat::dstep_red_mat(long im, long ido)
05935 {
05936 long i;
05937 double ret = 1.0;
05938 double dstep_r;
05939
05940 for (i=0; i<tnip; i++)
05941 {
05942 dstep_r = dstep_red_ip(i, 0, 0);
05943 if (dstep_r < ret)
05944 ret = dstep_r;
05945 }
05946 return ret;
05947 }
05948
05949
05950
05951
05952
05953
05954
05955
05956
05957
05958
05959
05960
05961
05962
05963
05964 double mechmat::dstep_red_ip(long ipp, long im, long ido)
05965 {
05966 double ret = 1.0;
05967
05968 switch (ip[ipp].tm[im])
05969 {
05970 case elisomat :
05971 case elortomat :
05972 case simplas1d:
05973 case jflow:
05974 case jflow2:
05975 case mohrcoul:
05976 case mohrcoulparab:
05977 case boermaterial:
05978 case druckerprager:
05979 case doubledrprager:
05980 case druckerprager2:
05981 case chenplast:
05982 case modcamclaymat:
05983 case modcamclaycoupmat:
05984 case hissplasmat:
05985 case microplaneM4:
05986 case microsimp:
05987 case microfibro:
05988 case shefpl:
05989 case glasgowdamage:
05990 case creep_damage:
05991 case shrinkagemat:
05992 case time_switchmat:
05993 case graphm:
05994 case varelisomat :
05995 case geoelast:
05996 case scaldamage:
05997 case fixortodamage:
05998 case scaldamagecc:
05999 case ortodamage:
06000 case ortodamagerot:
06001 case anisodamage:
06002 case anisodamagerot:
06003 case damage_plasticity:
06004 case simvisc:
06005 case lemaitr:
06006 case viscoplasticity:
06007 case creeprs:
06008 case creepb3:
06009 case creepdpl:
06010 case creepbaz:
06011 case creepeffym:
06012 case lenjonespot:
06013 case contmat:
06014 case cebfipcontmat:
06015 case consolidation:
06016 ret = 1.0;
06017 break;
06018
06019 case effective_stress:
06020 case nonlocdamgmat:
06021 case nonlocplastmat:
06022 case nonlocalmicroM4:
06023 ret = dstep_red_ip(ipp,im+1,ido);
06024 break;
06025
06026 case hypoplastmat:
06027 ret = hypopl[ip[ipp].idm[im]].dstep_red(ipp,im,ido);
06028 break;
06029 default:{
06030 print_err ("unknown material type is required",__FILE__,__LINE__,__func__);
06031 break;
06032 }
06033 }
06034 return ret;
06035 }
06036
06037
06038
06039
06040
06041
06042
06043
06044
06045
06046
06047
06048
06049
06050
06051
06052
06053
06054
06055
06056
06057
06058
06059
06060
06061
06062 double mechmat::yieldfunction (long ipp, long idpm, matrix &sig,vector &q)
06063 {
06064 double f;
06065
06066 switch (ip[ipp].tm[idpm]){
06067 case simplas1d:{
06068 f = spl1d[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06069 break;
06070 }
06071 case jflow:{
06072 f = j2f[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06073 break;
06074 }
06075 case jflow2:{
06076 f = j2f2[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06077 break;
06078 }
06079 case mohrcoulparab:{
06080 f = mcpar[ip[ipp].idm[idpm]].yieldfunction (sig);
06081 break;
06082 }
06083 case boermaterial:{
06084 f = boerm[ip[ipp].idm[idpm]].yieldfunction (sig);
06085 break;
06086 }
06087 case druckerprager:{
06088 f = drprm[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06089 break;
06090 }
06091 case doubledrprager:{
06092 f = ddpm[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06093 break;
06094 }
06095 case chenplast:{
06096 f = chplast[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06097 break;
06098 }
06099 case modcamclaymat:{
06100 f = cclay[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06101 break;
06102 }
06103 case modcamclaycoupmat:{
06104 f = cclayc[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06105 break;
06106 }
06107 case hissplasmat:{
06108 f = hisspl[ip[ipp].idm[idpm]].yieldfunction (sig,q);
06109 break;
06110 }
06111 default:
06112 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06113
06114 }
06115 return f;
06116 }
06117
06118
06119
06120
06121
06122
06123
06124
06125
06126
06127
06128
06129
06130
06131
06132
06133
06134
06135
06136 void mechmat::dfdsigma (long ipp, long idpm, matrix &sig, vector &q, vector &dfds)
06137 {
06138 matrix dfdst(3,3);
06139
06140 switch (ip[ipp].tm[idpm]){
06141 case simplas1d:{
06142 spl1d[ip[ipp].idm[idpm]].dfdsigma (sig,dfdst);
06143 break;
06144 }
06145 case jflow:{
06146 j2f[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
06147 break;
06148 }
06149 case jflow2:{
06150 j2f2[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
06151 break;
06152 }
06153 case mohrcoulparab:{
06154 mcpar[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
06155 break;
06156 }
06157 case boermaterial:{
06158 boerm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst);
06159 break;
06160 }
06161 case druckerprager:{
06162 drprm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst,q);
06163 break;
06164 }
06165 case doubledrprager:{
06166 ddpm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dfdst,q);
06167 break;
06168 }
06169 case chenplast:{
06170 chplast[ip[ipp].idm[idpm]].deryieldfdsigma (sig,q,dfdst);
06171 break;
06172 }
06173 case modcamclaymat:{
06174 cclay[ip[ipp].idm[idpm]].deryieldfsigma (sig, q, dfdst);
06175 break;
06176 }
06177 case modcamclaycoupmat:{
06178 cclayc[ip[ipp].idm[idpm]].deryieldfsigma (sig, q, dfdst);
06179 break;
06180 }
06181 default:
06182 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06183
06184 }
06185
06186
06187 tensor_vector (dfds,dfdst,ip[ipp].ssst,stress);
06188
06189 }
06190
06191
06192
06193
06194
06195
06196
06197
06198
06199
06200
06201
06202
06203
06204
06205
06206
06207
06208
06209 void mechmat::dgdsigma (long ipp, long idpm, matrix &sig, vector &q, vector &dgds)
06210 {
06211 matrix dgdst(3,3);
06212
06213 switch (ip[ipp].tm[idpm]){
06214 case simplas1d:{
06215 spl1d[ip[ipp].idm[idpm]].dfdsigma (sig,dgdst);
06216 break;
06217 }
06218 case jflow:{
06219 j2f[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst);
06220 break;
06221 }
06222 case jflow2:{
06223 j2f2[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst);
06224 break;
06225 }
06226 case mohrcoulparab:{
06227 mcpar[ip[ipp].idm[idpm]].derplaspotsigma (sig,dgdst);
06228 break;
06229 }
06230 case boermaterial:{
06231 boerm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst);
06232 break;
06233 }
06234 case druckerprager:{
06235 drprm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst,q);
06236 break;
06237 }
06238 case doubledrprager:{
06239 ddpm[ip[ipp].idm[idpm]].deryieldfsigma (sig,dgdst,q);
06240 break;
06241 }
06242 case chenplast:{
06243 chplast[ip[ipp].idm[idpm]].deryieldfdsigma (sig,q,dgdst);
06244 break;
06245 }
06246 case modcamclaymat:{
06247 cclay[ip[ipp].idm[idpm]].deryieldfsigma (sig, q, dgdst);
06248 break;
06249 }
06250 case modcamclaycoupmat:{
06251 cclayc[ip[ipp].idm[idpm]].deryieldfsigma (sig, q, dgdst);
06252 break;
06253 }
06254 default:
06255 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06256
06257 }
06258
06259
06260 tensor_vector (dgds,dgdst,ip[ipp].ssst,stress);
06261 }
06262
06263
06264
06265
06266
06267
06268
06269
06270
06271
06272
06273
06274
06275
06276
06277
06278
06279 void mechmat::dfdsigmadsigma (long ipp, long idpm, matrix &sig, matrix &dfdsds)
06280 {
06281 matrix dfdsdst(6,6);
06282
06283 switch (ip[ipp].tm[idpm]){
06284 case simplas1d:{
06285 spl1d[ip[ipp].idm[idpm]].dfdsigmadsigma (dfdsdst);
06286 break;
06287 }
06288 case jflow2:{
06289 j2f2[ip[ipp].idm[idpm]].deryielddsigmadsigma (sig,dfdsdst);
06290 break;
06291 }
06292 case modcamclaymat:{
06293 cclay[ip[ipp].idm[idpm]].dderyieldfsigma (dfdsdst);
06294 break;
06295 }
06296 case modcamclaycoupmat:{
06297 cclayc[ip[ipp].idm[idpm]].dderyieldfsigma (dfdsdst);
06298 break;
06299 }
06300 case chenplast:{
06301 chplast[ip[ipp].idm[idpm]].deryieldfdsigmadsigma (sig,dfdsdst);
06302 break;
06303 }
06304 default:
06305 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06306
06307 }
06308
06309
06310 tensor4_matrix (dfdsds,dfdsdst,ip[ipp].ssst);
06311
06312 }
06313
06314
06315
06316
06317
06318
06319
06320
06321
06322
06323
06324
06325
06326
06327
06328
06329
06330 void mechmat::dgdsigmadsigma (long ipp, long idpm, matrix &sig, matrix &dgdsds)
06331 {
06332 matrix dgdsdst(6,6);
06333
06334 switch (ip[ipp].tm[idpm]){
06335 case simplas1d:{
06336 spl1d[ip[ipp].idm[idpm]].dfdsigmadsigma (dgdsdst);
06337 break;
06338 }
06339 case jflow2:{
06340 j2f2[ip[ipp].idm[idpm]].deryielddsigmadsigma (sig,dgdsdst);
06341 break;
06342 }
06343 case modcamclaymat:{
06344 cclay[ip[ipp].idm[idpm]].dderyieldfsigma (dgdsdst);
06345 break;
06346 }
06347 case modcamclaycoupmat:{
06348 cclayc[ip[ipp].idm[idpm]].dderyieldfsigma (dgdsdst);
06349 break;
06350 }
06351 case chenplast:{
06352 chplast[ip[ipp].idm[idpm]].deryieldfdsigmadsigma (sig,dgdsdst);
06353 break;
06354 }
06355 default:
06356 print_err("unknown plasticity model is required", __FILE__, __LINE__, __func__);
06357
06358 }
06359
06360
06361 tensor4_matrix (dgdsds,dgdsdst,ip[ipp].ssst);
06362
06363 }
06364
06365
06366
06367
06368
06369
06370
06371
06372
06373
06374
06375
06376
06377
06378
06379
06380
06381
06382 void mechmat::dfdqpar (long ipp, long idpm,matrix &sig,vector &q,vector &dq)
06383 {
06384 switch (ip[ipp].tm[idpm]){
06385 case simplas1d:{
06386 spl1d[ip[ipp].idm[idpm]].dfdqpar (dq);
06387 break;
06388 }
06389 case jflow:{
06390 j2f[ip[ipp].idm[idpm]].deryieldfq (dq);
06391 break;
06392 }
06393 case jflow2:{
06394 j2f2[ip[ipp].idm[idpm]].deryieldfq (dq);
06395 break;
06396 }
06397 case mohrcoulparab:{
06398 break;
06399 }
06400 case boermaterial:{
06401 break;
06402 }
06403 case druckerprager:{
06404 break;
06405 }
06406 case doubledrprager:{
06407 break;
06408 }
06409 case modcamclaymat:{
06410 cclay[ip[ipp].idm[idpm]].deryieldfq (sig, dq);
06411 break;
06412 }
06413 case modcamclaycoupmat:{
06414 cclayc[ip[ipp].idm[idpm]].deryieldfq (sig, q, dq);
06415 break;
06416 }
06417 case chenplast:{
06418 chplast[ip[ipp].idm[idpm]].deryieldfdq (sig,q,dq);
06419 break;
06420 }
06421 default:
06422 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06423
06424 }
06425 }
06426
06427
06428
06429
06430
06431
06432
06433
06434
06435
06436
06437
06438
06439
06440
06441
06442
06443
06444 void mechmat::dgdqpar (long ipp, long idpm,matrix &sig,vector &q,vector &dq)
06445 {
06446 switch (ip[ipp].tm[idpm]){
06447 case simplas1d:{
06448 spl1d[ip[ipp].idm[idpm]].dfdqpar (dq);
06449 break;
06450 }
06451 case jflow:{
06452 j2f[ip[ipp].idm[idpm]].deryieldfq (dq);
06453 break;
06454 }
06455 case jflow2:{
06456 j2f2[ip[ipp].idm[idpm]].deryieldfq (dq);
06457 break;
06458 }
06459 case mohrcoulparab:{
06460 break;
06461 }
06462 case boermaterial:{
06463 break;
06464 }
06465 case druckerprager:{
06466 break;
06467 }
06468 case doubledrprager:{
06469 break;
06470 }
06471 case modcamclaymat:{
06472 cclay[ip[ipp].idm[idpm]].deryieldfq (sig, dq);
06473 break;
06474 }
06475 case modcamclaycoupmat:{
06476 cclayc[ip[ipp].idm[idpm]].deryieldfq (sig, q, dq);
06477 break;
06478 }
06479 case chenplast:{
06480 chplast[ip[ipp].idm[idpm]].deryieldfdq (sig,q,dq);
06481 break;
06482 }
06483 default:
06484 print_err ("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06485
06486 }
06487 }
06488
06489
06490
06491
06492
06493
06494
06495
06496
06497
06498
06499
06500
06501
06502
06503
06504
06505
06506
06507
06508 void mechmat::dfdsigmadq (long ipp, long idpm,matrix &sig,vector &q,matrix &dfdsdq)
06509 {
06510
06511 long nh=dfdsdq.n;
06512
06513 matrix dfdsdqt(6,nh);
06514
06515 switch (ip[ipp].tm[idpm]){
06516 case simplas1d:{
06517 spl1d[ip[ipp].idm[idpm]].dfdsigmadq (dfdsdqt);
06518 break;
06519 }
06520 case jflow2:{
06521 j2f2[ip[ipp].idm[idpm]].deryieldfdsigmadq (dfdsdqt);
06522 break;
06523 }
06524 case modcamclaymat:{
06525 cclay[ip[ipp].idm[idpm]].deryieldfdsdq (dfdsdqt);
06526 break;
06527 }
06528 case modcamclaycoupmat:{
06529 cclayc[ip[ipp].idm[idpm]].deryieldfdsdq (dfdsdqt);
06530 break;
06531 }
06532 case chenplast:{
06533 chplast[ip[ipp].idm[idpm]].deryieldfdsigmadq (sig,q,dfdsdqt);
06534 break;
06535 }
06536 default:
06537 print_err("unknown plasticity model is required", __FILE__,__LINE__,__func__);
06538
06539 }
06540
06541
06542 gentensor_matrix (dfdsdq,dfdsdqt,ip[ipp].ssst);
06543
06544 }
06545
06546
06547
06548
06549
06550
06551
06552
06553
06554
06555
06556
06557
06558
06559
06560
06561
06562
06563
06564
06565 void mechmat::dgdsigmadq (long ipp, long idpm,matrix &sig,vector &q,matrix &dgdsdq)
06566 {
06567
06568 long nh=dgdsdq.n;
06569
06570 matrix dgdsdqt(6,nh);
06571
06572 switch (ip[ipp].tm[idpm]){
06573 case simplas1d:{
06574 spl1d[ip[ipp].idm[idpm]].dfdsigmadq (dgdsdqt);
06575 break;
06576 }
06577 case jflow2:{
06578 j2f2[ip[ipp].idm[idpm]].deryieldfdsigmadq (dgdsdqt);
06579 break;
06580 }
06581 case modcamclaymat:{
06582 cclay[ip[ipp].idm[idpm]].deryieldfdsdq (dgdsdqt);
06583 break;
06584 }
06585 case modcamclaycoupmat:{
06586 cclayc[ip[ipp].idm[idpm]].deryieldfdsdq (dgdsdqt);
06587 break;
06588 }
06589 case chenplast:{
06590 chplast[ip[ipp].idm[idpm]].deryieldfdsigmadq (sig,q,dgdsdqt);
06591 break;
06592 }
06593 default:
06594 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06595
06596 }
06597
06598
06599 gentensor_matrix (dgdsdq,dgdsdqt,ip[ipp].ssst);
06600 }
06601
06602
06603
06604
06605
06606
06607
06608
06609
06610
06611
06612
06613
06614
06615
06616
06617 void mechmat::dhdsigma (long ipp,long idpm,matrix &sigt,vector &q,matrix &dhds)
06618 {
06619
06620 long nh=dhds.m;
06621
06622 matrix dhdst(6,nh);
06623
06624 switch (ip[ipp].tm[idpm]){
06625 case simplas1d:{
06626
06627 break;
06628 }
06629 case chenplast:{
06630 chplast[ip[ipp].idm[idpm]].dhdsigma (sigt,q,dhdst);
06631 break;
06632 }
06633 default:
06634 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06635
06636 }
06637
06638
06639 gentensor_matrix (dhds,dhdst,ip[ipp].ssst);
06640 }
06641
06642
06643
06644
06645
06646
06647
06648
06649
06650
06651
06652
06653
06654
06655 void mechmat::dhdgamma (long ipp,long idpm,vector &epsp,matrix &sig,vector &q,vector &dhdc)
06656 {
06657 switch (ip[ipp].tm[idpm]){
06658 case simplas1d:{
06659 spl1d[ip[ipp].idm[idpm]].dhdgamma (ipp,epsp,sig,dhdc);
06660 break;
06661 }
06662 case chenplast:{
06663 chplast[ip[ipp].idm[idpm]].dhdgamma (dhdc);
06664 break;
06665 }
06666 default:
06667 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06668
06669 }
06670 }
06671
06672
06673
06674
06675
06676
06677
06678
06679
06680
06681
06682
06683
06684
06685
06686
06687
06688
06689 void mechmat::dhdqpar (long ipp,long idpm,long ido,matrix &sigt,double dgamma, vector &q,matrix &dhdq)
06690 {
06691 switch (ip[ipp].tm[idpm]){
06692 case simplas1d:{
06693
06694 break;
06695 }
06696 case chenplast:{
06697 chplast[ip[ipp].idm[idpm]].dhdqpar (sigt,q,dhdq);
06698 break;
06699 }
06700 case modcamclaymat:{
06701 cclay[ip[ipp].idm[idpm]].dhardfdq (ipp, ido, dgamma, q, dhdq);
06702 break;
06703 }
06704 default:
06705 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06706
06707 }
06708 }
06709
06710
06711
06712
06713
06714
06715
06716
06717
06718
06719
06720
06721
06722
06723
06724
06725 void mechmat::hardsoftfunction (long ipp,long idpm,matrix &sigt,vector &q,vector &h)
06726 {
06727 switch (ip[ipp].tm[idpm]){
06728 case chenplast:{
06729 chplast[ip[ipp].idm[idpm]].hvalues (sigt,q,h);
06730 break;
06731 }
06732 default:
06733 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06734
06735 }
06736 }
06737
06738
06739
06740
06741
06742
06743
06744
06745
06746
06747
06748
06749
06750
06751
06752
06753
06754 void mechmat::plasticmoduli (long ipp,long idpm,vector &epsp,matrix &sig,matrix &h)
06755 {
06756 switch (ip[ipp].tm[idpm]){
06757 case simplas1d:{
06758 spl1d[ip[ipp].idm[idpm]].plasmod (h);
06759 break;
06760 }
06761 case jflow:{
06762 j2f[ip[ipp].idm[idpm]].plasmod (h);
06763 break;
06764 }
06765 case jflow2:{
06766 j2f2[ip[ipp].idm[idpm]].plasmod (h);
06767 break;
06768 }
06769 case chenplast:{
06770
06771 break;
06772 }
06773 default:
06774 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06775
06776 }
06777
06778 }
06779
06780
06781
06782
06783
06784
06785
06786
06787
06788
06789
06790
06791
06792
06793
06794
06795
06796
06797
06798
06799
06800
06801
06802
06803
06804
06805
06806
06807
06808
06809
06810
06811
06812
06813
06814
06815
06816
06817
06818
06819
06820
06821
06822
06823
06824 void mechmat::plasmodscalar(double &denom,long ipp, long idpm, long ido, matrix &sig, vector &eps, vector &epsp, vector &qtr,double gamma)
06825 {
06826 double ret = 0.0;
06827 switch (ip[ipp].tm[idpm]){
06828 case simplas1d:{
06829
06830 ret = spl1d[ip[ipp].idm[idpm]].plasmodscalar (sig,epsp,qtr,gamma);
06831 break;
06832 }
06833 case jflow:{
06834 ret = j2f[ip[ipp].idm[idpm]].plasmodscalar (qtr);
06835 break;
06836 }
06837 case jflow2:{
06838 ret = j2f2[ip[ipp].idm[idpm]].plasmodscalar (qtr);
06839 break;
06840 }
06841 case mohrcoulparab:{
06842 break;
06843 }
06844 case boermaterial:{
06845 break;
06846 }
06847 case druckerprager:{
06848 ret = drprm[ip[ipp].idm[idpm]].plasmodscalar(qtr);
06849 break;
06850 }
06851 case doubledrprager:{
06852 ret = ddpm[ip[ipp].idm[idpm]].plasmodscalar(sig, qtr);
06853 break;
06854 }
06855 case chenplast:{
06856
06857 if (ret<1.0e-6)
06858 denom*=2.0;
06859
06860
06861
06862
06863
06864 return;
06865
06866 }
06867 case modcamclaymat:{
06868 ret = cclay[ip[ipp].idm[idpm]].plasmodscalar(ipp, ido, sig, qtr);
06869 break;
06870 }
06871 case modcamclaycoupmat:{
06872 ret = cclayc[ip[ipp].idm[idpm]].plasmodscalar(ipp, ido, sig, qtr);
06873 break;
06874 }
06875 default:
06876 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06877
06878 }
06879 denom = ret;
06880
06881 }
06882
06883
06884
06885
06886
06887
06888
06889
06890
06891
06892
06893
06894
06895 long mechmat::give_num_interparam (long ipp,long idpm)
06896 {
06897 long np;
06898
06899 switch (ip[ipp].tm[idpm]){
06900 case simplas1d:{
06901 np = spl1d[ip[ipp].idm[idpm]].give_num_interparam ();
06902 break;
06903 }
06904 case jflow:{
06905 np = j2f[ip[ipp].idm[idpm]].give_num_interparam ();
06906 break;
06907 }
06908 case mohrcoul:{
06909 break;
06910 }
06911 case mohrcoulparab:{
06912 break;
06913 }
06914 case boermaterial:{
06915 break;
06916 }
06917 case druckerprager:{
06918 break;
06919 }
06920 case modcamclaymat:{
06921 break;
06922 }
06923 case modcamclaycoupmat:{
06924 break;
06925 }
06926 default:
06927 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
06928
06929 }
06930
06931 return np;
06932 }
06933
06934
06935
06936
06937
06938
06939
06940
06941
06942
06943
06944
06945
06946
06947
06948 void mechmat::give_interparam (long ipp,long im,long ido,vector &q)
06949 {
06950 switch (ip[ipp].tm[im]){
06951 case simplas1d:{
06952 spl1d[ip[ipp].idm[im]].give_interparam (ipp,ido,q);
06953 break;
06954 }
06955 case jflow:{
06956 j2f[ip[ipp].idm[im]].give_interparam (ipp,ido,q);
06957 break;
06958 }
06959 case mohrcoul:{
06960 break;
06961 }
06962 case mohrcoulparab:{
06963 break;
06964 }
06965 case boermaterial:{
06966 break;
06967 }
06968 case druckerprager:{
06969 break;
06970 }
06971 case modcamclaymat:{
06972 break;
06973 }
06974 case modcamclaycoupmat:{
06975 break;
06976 }
06977
06978
06979
06980
06981
06982 default:
06983 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
06984
06985 }
06986
06987 }
06988
06989
06990
06991
06992
06993
06994
06995
06996
06997
06998
06999
07000
07001
07002
07003
07004
07005
07006
07007
07008 void mechmat::updateq (long ipp, long idpm, long ido,double dgamma, vector &eps, vector &epsp, matrix &sig, vector &q)
07009 {
07010 switch (ip[ipp].tm[idpm]){
07011 case simplas1d:{
07012
07013 spl1d[ip[ipp].idm[idpm]].updateq (ipp,dgamma, epsp,q);
07014 break;
07015 }
07016 case jflow:{
07017 j2f[ip[ipp].idm[idpm]].updateq (dgamma, q);
07018 break;
07019 }
07020 case jflow2:{
07021 j2f2[ip[ipp].idm[idpm]].updateq (dgamma, q);
07022 break;
07023 }
07024 case mohrcoul:{
07025 break;
07026 }
07027 case mohrcoulparab:{
07028 break;
07029 }
07030 case boermaterial:{
07031 break;
07032 }
07033 case druckerprager:{
07034 drprm[ip[ipp].idm[idpm]].updateq(ipp, epsp, q);
07035 break;
07036 }
07037 case doubledrprager:{
07038 ddpm[ip[ipp].idm[idpm]].updateq(ipp, dgamma, q);
07039 break;
07040 }
07041 case chenplast:{
07042
07043 break;
07044 }
07045 case modcamclaymat:{
07046 cclay[ip[ipp].idm[idpm]].updateq(ipp, ido, eps, epsp, sig, q);
07047 break;
07048 }
07049 case modcamclaycoupmat:{
07050 cclayc[ip[ipp].idm[idpm]].updateq(ipp, ido, eps, epsp, sig, q);
07051 break;
07052 }
07053 default:
07054 print_err("unknown plasticity model is required",__FILE__,__LINE__,__func__);
07055
07056 }
07057 }
07058
07059
07060
07061
07062
07063
07064
07065
07066
07067
07068
07069
07070
07071
07072
07073
07074
07075
07076
07077
07078
07079
07080
07081
07082
07083
07084 void mechmat::cutting_plane (long ipp,long im,long ido,double &gamma,vector &epsn,vector &epsp,vector &q,long ni,double err)
07085 {
07086 long i,j,ncomp=epsn.n,idem;
07087 double f,denom,denomh,dgamma,nu;
07088 vector epsa(ncomp),sig(ncomp),dfds(ncomp),dgds(ncomp);
07089 matrix d(ncomp,ncomp),sigt(3,3);
07090
07091
07092 dgamma=0.0;
07093
07094
07095 elmatstiff (d,ipp);
07096 if (ip[ipp].ssst == planestrain)
07097 {
07098 d[0][3] = d[0][1]; d[1][3] = d[1][0];
07099 d[3][0] = d[1][0]; d[3][1] = d[1][0]; d[3][3] = d[1][1];
07100 }
07101 idem = ip[ipp].gemid();
07102 if (ip[ipp].ssst == planestress)
07103 {
07104 nu = Mm->give_actual_nu(ipp);
07105 epsn[3] = -nu / (1.0 - nu) * (epsn[0]-epsp[0]+epsn[1]-epsp[1]);
07106 }
07107
07108
07109
07110
07111
07112 for (i=0;i<ni;i++){
07113
07114
07115
07116
07117 subv (epsn,epsp,epsa);
07118
07119 mxv (d,epsa,sig);
07120
07121 vector_tensor (sig,sigt,ip[ipp].ssst,stress);
07122
07123 updateq(ipp, im, ido, dgamma, epsn, epsp, sigt, q);
07124
07125 f = yieldfunction (ipp,im,sigt,q);
07126
07127 if (f<err) break;
07128 if (i==ni-1 && f>err){
07129 print_err("yield surface was not reached",__FILE__,__LINE__,__func__);
07130 fprintf (Out,"\n yield surface was not reached");
07131 }
07132
07133 dfdsigma (ipp, im, sigt, q, dfds);
07134 dgdsigma (ipp, im, sigt, q, dgds);
07135
07136 mxv (d,dgds,epsa);
07137 scprd (dfds,epsa,denom);
07138
07139 denomh=0;
07140 plasmodscalar(denomh,ipp, im, ido, sigt, epsn, epsp, q,gamma);
07141 denom += denomh;
07142
07143
07144 dgamma = f/denom;
07145
07146 gamma+=dgamma;
07147
07148 for (j=0;j<ncomp;j++){
07149 epsp[j]+=dgamma*dgds[j];
07150 }
07151 if (ip[ipp].ssst == planestress)
07152 epsp[3] = epsp[0]+epsp[1];
07153 }
07154
07155
07156 subv (epsn,epsp,epsa);
07157
07158 mxv (d,epsa,sig);
07159
07160 storestress (0,ipp,sig);
07161 }
07162
07163
07164
07165
07166
07167
07168
07169
07170
07171
07172
07173
07174
07175
07176
07177
07178
07179
07180
07181
07182
07183
07184
07185
07186
07187
07188
07189 void mechmat::cutting_plane2 (long ipp,long im,long ido,double &gamma,vector &epsn,vector &epsp,vector &q,long ni,double err)
07190 {
07191 long i,j,b,ncomp=epsn.n,idem;
07192 double f,denom,denomh,dgamma,pom;
07193 vector epsa(ncomp),sig(ncomp),dfds(ncomp),dgds(ncomp);
07194 matrix d(ncomp,ncomp),sigt(3,3);
07195
07196
07197 dgamma=0.0;
07198
07199
07200 if (ip[ipp].ssst == planestress)
07201 {
07202 ip[ipp].ssst = axisymm;
07203 b=1;
07204 }
07205
07206
07207 elmatstiff (d,ipp);
07208
07209 if (ip[ipp].ssst == planestrain)
07210 {
07211 d[0][3] = d[0][1]; d[1][3] = d[1][0];
07212 d[3][0] = d[1][0]; d[3][1] = d[1][0]; d[3][3] = d[1][1];
07213 }
07214 idem = ip[ipp].gemid();
07215
07216
07217
07218
07219 for (i=0;i<ni;i++){
07220
07221
07222
07223
07224 subv (epsn,epsp,epsa);
07225
07226 mxv (d,epsa,sig);
07227
07228
07229
07230
07231 vector_tensor (sig,sigt,ip[ipp].ssst,stress);
07232
07233 updateq(ipp, im, ido, dgamma, epsn, epsp, sigt, q);
07234
07235 f = yieldfunction (ipp,im,sigt,q);
07236
07237
07238
07239
07240 if (f<err) break;
07241 if (i==ni-1 && f>err){
07242 print_err("yield surface was not reached",__FILE__,__LINE__,__func__);
07243 fprintf (Out,"\n yield surface was not reached");
07244 }
07245
07246 dfdsigma (ipp, im, sigt, q, dfds);
07247 dgdsigma (ipp, im, sigt, q, dgds);
07248
07249 mxv (d,dgds,epsa);
07250 scprd (dfds,epsa,denom);
07251
07252
07253
07254
07255 denomh=0;
07256 plasmodscalar(denomh,ipp, im, ido, sigt, epsn, epsp, q,gamma);
07257 denom += denomh;
07258
07259
07260
07261
07262 dgamma = f/denom;
07263
07264 gamma+=dgamma;
07265
07266 for (j=0;j<ncomp;j++){
07267 epsp[j]+=dgamma*dgds[j];
07268 }
07269 }
07270
07271
07272 subv (epsn,epsp,epsa);
07273
07274 mxv (d,epsa,sig);
07275 if (b==1)
07276 {
07277 pom = sig[3];
07278 sig[3] = sig[2];
07279 sig[2] = pom;
07280 }
07281
07282 storestress (0,ipp,sig);
07283
07284 if (b==1)
07285 ip[ipp].ssst = planestress;
07286 }
07287
07288
07289
07290
07291
07292
07293
07294
07295
07296
07297
07298
07299
07300
07301
07302
07303
07304
07305
07306
07307
07308
07309
07310
07311
07312 void mechmat::mult_surf_cutting_plane (long ipp,vector &epsn,vector &epsp,vector &gamma,vector &q)
07313 {
07314 long ni,i,j,k,nas,ncomp=epsn.n,nhp=q.n,nyf=gamma.n,*stat;
07315 double err;
07316 vector epsptr(ncomp),epsa(ncomp),sig(ncomp),f(nyf),qtr(nhp);
07317 vector dq,dgamma,ff;
07318 matrix d(ncomp,ncomp),h(nhp,nhp);
07319 matrix fsig,fq,cpm,amq,am,acpm;
07320
07321
07322
07323 ni=0;
07324 err=10.0;
07325
07326
07327 copyv (epsp,epsptr);
07328 copyv (q,qtr);
07329
07330 stat = new long [nyf];
07331
07332
07333 for (i=0;i<ni;i++){
07334
07335 subv (epsn,epsptr,epsa);
07336
07337
07338 elmatstiff (d,ipp);
07339
07340
07341 mxv (d,epsa,sig);
07342
07343
07344
07345
07346
07347
07348
07349
07350
07351
07352
07353
07354
07355
07356
07357
07358 nas=0;
07359 for (j=0;j<nyf;j++){
07360 if (f[j]<err){ stat[j]=0; }
07361 else{ nas++; stat[j]=1; }
07362 }
07363 if (nas==0) break;
07364
07365
07366 allocm (nas,6,fsig); allocm (nas,6,am); allocm (nas,nas,cpm);
07367 allocm (nas,nhp,fq); allocm (nhp,nas,amq); allocm (nas,nas,acpm);
07368 allocv (nas,ff); allocv (nas,dgamma); allocv (nas,dq);
07369
07370
07371
07372
07373
07374
07375
07376
07377
07378
07379
07380
07381
07382
07383
07384
07385
07386
07387
07388
07389
07390
07391 mxm (fsig,d,am);
07392 mxmt (am,fsig,cpm);
07393
07394
07395 mxmt (h,fq,amq);
07396 mxm (fq,amq,acpm);
07397
07398 addm (cpm,acpm,cpm);
07399
07400
07401 gemp (cpm.a,dgamma.a,ff.a,nas,1,Mp->zero,1);
07402
07403
07404 k=0;
07405 for (j=0;j<nyf;j++){
07406 if (stat[j]==1){
07407 gamma[j]+=dgamma[k];
07408 k++;
07409 }
07410 }
07411
07412 mtxv (fsig,dgamma,sig);
07413 mxv (amq,dgamma,dq);
07414
07415
07416 for (j=0;j<ncomp;j++){
07417 epsptr[j]+=sig[j];
07418 }
07419
07420 for (j=0;j<nhp;j++){
07421 qtr[j]-=dq[j];
07422 }
07423
07424
07425 destrm (fsig); destrm (am); destrm (cpm);
07426 destrm (fq); destrm (amq); destrm (acpm);
07427 destrv (ff); destrv (dgamma); destrv (dq);
07428
07429 }
07430
07431 copyv (epsptr,epsp);
07432 copyv (qtr,q);
07433
07434
07435
07436 subv (epsn,epsp,epsa);
07437
07438
07439 elmatstiff (d,ipp);
07440
07441
07442 mxv (d,epsa,sig);
07443
07444 storestress (0,ipp,sig);
07445
07446 }
07447
07448
07449
07450
07451
07452
07453
07454
07455
07456
07457
07458
07459
07460
07461
07462
07463
07464
07465
07466
07467
07468
07469
07470
07471
07472 void mechmat::newton_stress_return (long ipp,long im,long ido,double &gamma,vector &epsn,vector &epsp,vector &q,long ni,double err)
07473 {
07474 long i,j,k,n,ncompstr,ncomphard;
07475 double zero,f,dgamma,norstr,norstrincr,norconspar,norhardpar,norhardparincr;
07476
07477
07478 zero=Mp->zero;
07479
07480
07481 ncompstr = epsn.n;
07482
07483 ncomphard = q.n;
07484
07485 n = ncompstr + 1 + ncomphard;
07486
07487
07488
07489
07490
07491
07492
07493
07494 double *a;
07495
07496 double *x;
07497
07498 double *y;
07499
07500 matrix d(ncompstr,ncompstr);
07501
07502 matrix c(ncompstr,ncompstr);
07503
07504 vector epse(ncompstr);
07505
07506 vector sig(ncompstr);
07507
07508 matrix sigt(3,3);
07509
07510 vector sigtrial(ncompstr);
07511
07512
07513
07514 matrix dgdsds(ncompstr,ncompstr);
07515
07516 vector dgds(ncompstr);
07517
07518 vector dfds(ncompstr);
07519
07520 matrix dgdsdq(ncompstr,ncomphard);
07521
07522 vector dfdq(ncomphard);
07523
07524 matrix dhds (ncomphard,ncompstr);
07525
07526 vector dhdc (ncomphard);
07527
07528 matrix dhdq(ncomphard,ncomphard);
07529
07530 vector hs(ncomphard);
07531
07532
07533 matrix auxm(ncompstr,ncompstr),auxm2(ncompstr,ncomphard);
07534
07535 vector auxv(ncompstr),auxv2(ncompstr);
07536
07537
07538
07539 a = new double [n*n];
07540 x = new double [n];
07541 y = new double [n];
07542
07543
07544 elmatstiff (d,ipp);
07545
07546 elmatcompl (c,ipp);
07547
07548
07549 dgamma=0.0;
07550
07551
07552 subv (epsn,epsp,epse);
07553
07554 mxv (d,epse,sig);
07555
07556
07557 if (ic){
07558 for (j=0; j < ncompstr; j++)
07559 sig[j] += ic[ipp][3+j];
07560 }
07561
07562
07563 vector_tensor (sig,sigt,ip[ipp].ssst,stress);
07564
07565
07566
07567
07568 f = yieldfunction (ipp,im,sigt,q);
07569
07570
07571 hardsoftfunction (ipp,im,sigt,q,hs);
07572
07573 if (f>err){
07574
07575
07576
07577
07578 copyv (sig,sigtrial);
07579
07580
07581 for (i=0;i<ni;i++){
07582
07583
07584
07585 nullv (a,n*n);
07586 nullv (x,n);
07587 nullv (y,n);
07588
07589
07590
07591
07592
07593 dgdsigmadsigma (ipp, im, sigt, dgdsds);
07594
07595 mxm (d,dgdsds,auxm);
07596 cmulm (dgamma,auxm);
07597
07598
07599
07600 for (j=0;j<ncompstr;j++){
07601 for (k=0;k<ncompstr;k++){
07602 a[j*n+k]=auxm[j][k];
07603 }
07604 a[j*n+j]+=1.0;
07605 }
07606
07607
07608
07609
07610 dgdsigma (ipp, im, sigt, q, dgds);
07611
07612 mxv (d,dgds,auxv);
07613
07614
07615
07616 for (j=0;j<ncompstr;j++){
07617 a[j*n+ncompstr]=auxv[j];
07618 }
07619
07620
07621
07622
07623 dfdsigma (ipp, im, sigt, q, dfds);
07624
07625
07626
07627 for (j=0;j<ncompstr;j++){
07628 a[ncompstr*n+j]=dfds[j];
07629 }
07630
07631
07632
07633
07634
07635
07636 if (ncomphard>0){
07637
07638
07639
07640
07641
07642
07643
07644
07645 for (j=0;j<ncompstr;j++){
07646 for (k=0;k<ncomphard;k++){
07647 a[j*n+ncompstr+1+k]=auxm2[j][k]*dgamma;
07648 }
07649 }
07650
07651
07652
07653
07654
07655
07656
07657
07658
07659
07660 for (j=0;j<ncomphard;j++){
07661 for (k=0;k<ncompstr;k++){
07662 a[(ncompstr+1+j)*n+k]=dhds[j][k];
07663 }
07664 }
07665
07666
07667
07668
07669
07670
07671
07672
07673 for (j=0;j<ncomphard;j++){
07674 a[ncompstr*n+ncompstr+1+j]=dfdq[j];
07675 }
07676
07677
07678
07679
07680
07681
07682
07683
07684 for (j=0;j<ncomphard;j++){
07685 a[(ncompstr+1+j)*n+ncompstr]=0.0-dhdc[j];
07686 }
07687
07688
07689
07690
07691
07692
07693
07694
07695
07696 for (j=0;j<ncomphard;j++){
07697 for (k=0;k<ncomphard;k++){
07698 a[(ncompstr+1+j)*n+ncompstr+1+k]=dhdq[j][k];
07699 }
07700 a[(ncompstr+1+j)*n+ncompstr+1+j]-=1.0;
07701 }
07702 }
07703
07704
07705
07706
07707
07708 for (j=0;j<ncompstr;j++){
07709 y[j]=sigtrial[j]-sig[j]-dgamma*auxv[j];
07710 }
07711
07712
07713 y[ncompstr]=0.0-f;
07714
07715
07716 for (j=0;j<ncomphard;j++){
07717
07718
07719 }
07720
07721
07722
07723
07724
07725
07726
07727
07728
07729
07730
07731
07732
07733
07734
07735 gemp (a,x,y,n,1,1.0e-10,1);
07736
07737
07738
07739
07740
07741
07742
07743 norstr=0.0;
07744
07745 norstrincr=0.0;
07746
07747 for (j=0;j<ncompstr;j++){
07748 norstrincr+=x[j]*x[j];
07749 norstr+=sig[j]*sig[j];
07750 }
07751
07752 norconspar=x[ncompstr]*x[ncompstr];
07753
07754
07755 norhardpar=0.0;
07756
07757 norhardparincr=0.0;
07758
07759 for (j=0;j<ncomphard;j++){
07760 norhardpar+=q[j]*q[j];
07761 }
07762 for (j=ncompstr+1;j<ncompstr+1+ncomphard;j++){
07763 norhardparincr+=x[j]*x[j];
07764 }
07765
07766
07767
07768 for (j=0;j<ncompstr;j++){
07769 sig[j]+=x[j];
07770 auxv[j]=x[j];
07771 }
07772
07773
07774 dgamma+=x[ncompstr];
07775
07776
07777 for (j=0;j<ncomphard;j++){
07778 q[j]+=x[ncompstr+1+j];
07779 }
07780
07781
07782
07783
07784
07785
07786
07787
07788
07789
07790
07791
07792 mxv (c,auxv,auxv2);
07793
07794 for (j=0;j<ncompstr;j++){
07795 epsp[j]-=auxv2[j];
07796 }
07797
07798
07799
07800 vector_tensor (sig,sigt,ip[ipp].ssst,stress);
07801
07802
07803
07804
07805 f = yieldfunction (ipp,im,sigt,q);
07806
07807
07808 hardsoftfunction (ipp,im,sigt,q,hs);
07809
07810
07811
07812
07813
07814
07815
07816
07817
07818 j=0; k=0;
07819 if (norstr>zero){
07820 j++;
07821 if (norstrincr/norstr<err) k++;
07822 }
07823 if (norhardpar>zero){
07824 j++;
07825 if (norhardparincr/norhardpar<err) k++;
07826 }
07827 if (dgamma>zero){
07828 j++;
07829 if (norconspar/dgamma<err) k++;
07830 }
07831
07832
07833 if (j==k){
07834
07835
07836
07837
07838 break;
07839 }
07840 }
07841 }
07842
07843
07844 gamma+=dgamma;
07845
07846
07847 storestress (0,ipp,sig);
07848
07849 delete [] a;
07850 delete [] x;
07851 delete [] y;
07852 }
07853
07854
07855
07856
07857
07858
07859
07860
07861
07862
07863
07864
07865
07866
07867
07868
07869
07870
07871
07872
07873
07874
07875
07876
07877
07878 void mechmat::newton_stress_return_2 (long ipp,long im,long ido,double &gamma,vector &epsn,vector &epsp,vector &q,long ni,double err)
07879 {
07880 long i,j,k,n,ncompstr,ncomphard;
07881 double zero,f,dgamma,norstr,norstrincr,norconspar,norhardpar,norhardparincr;
07882
07883
07884
07885 zero=Mp->zero;
07886
07887
07888 ncompstr = epsn.n;
07889
07890 ncomphard = q.n;
07891
07892 n = ncompstr + 1 + ncomphard;
07893
07894
07895
07896
07897
07898
07899 vector qold(ncomphard);
07900
07901
07902 double *a;
07903
07904 double *x;
07905
07906 double *y;
07907
07908 matrix d(ncompstr,ncompstr);
07909
07910 matrix c(ncompstr,ncompstr);
07911
07912 vector epse(ncompstr);
07913
07914 vector sig(ncompstr);
07915
07916 matrix sigt(3,3);
07917
07918 vector sigtrial(ncompstr);
07919
07920
07921
07922 matrix dgdsds(ncompstr,ncompstr);
07923
07924 vector dgds(ncompstr);
07925
07926 vector dfds(ncompstr);
07927
07928 matrix dgdsdq(ncompstr,ncomphard);
07929
07930 vector dfdq(ncomphard);
07931
07932 matrix dhds (ncompstr,ncomphard);
07933
07934 vector dhdc (ncomphard);
07935
07936 matrix dhdq(ncomphard,ncomphard);
07937
07938 vector hs(ncomphard);
07939
07940
07941 matrix auxm(ncompstr,ncompstr),auxm2(ncompstr,ncomphard);
07942
07943 vector auxv(ncompstr),auxv2(ncompstr);
07944
07945
07946
07947 a = new double [n*n];
07948 x = new double [n];
07949 y = new double [n];
07950
07951
07952 elmatstiff (d,ipp);
07953
07954 elmatcompl (c,ipp);
07955
07956
07957 dgamma=0.0;
07958
07959
07960 subv (epsn,epsp,epse);
07961
07962 mxv (d,epse,sig);
07963
07964
07965 if (ic){
07966 for (j=0; j < ncompstr; j++)
07967 sig[j] += ic[ipp][3+j];
07968 }
07969
07970
07971 vector_tensor (sig,sigt,ip[ipp].ssst,stress);
07972
07973 copyv (q,qold);
07974
07975
07976 f = yieldfunction (ipp,im,sigt,q);
07977
07978
07979 hardsoftfunction (ipp,im,sigt,q,hs);
07980
07981 if (f>err){
07982
07983
07984
07985
07986 copyv (sig,sigtrial);
07987
07988
07989 for (i=0;i<ni;i++){
07990
07991
07992
07993 nullv (a,n*n);
07994 nullv (x,n);
07995 nullv (y,n);
07996
07997
07998
07999
08000
08001 dgdsigmadsigma (ipp, im, sigt, dgdsds);
08002
08003 mxm (d,dgdsds,auxm);
08004 cmulm (dgamma,auxm);
08005
08006
08007
08008 for (j=0;j<ncompstr;j++){
08009 for (k=0;k<ncompstr;k++){
08010 a[j*n+k]=auxm[j][k];
08011 }
08012 a[j*n+j]+=1.0;
08013 }
08014
08015
08016
08017
08018 dgdsigma (ipp, im, sigt, q, dgds);
08019
08020 mxv (d,dgds,auxv);
08021
08022
08023
08024 for (j=0;j<ncompstr;j++){
08025 a[j*n+ncompstr]=auxv[j];
08026 }
08027
08028
08029
08030
08031 dfdsigma (ipp, im, sigt, q, dfds);
08032
08033
08034
08035 for (j=0;j<ncompstr;j++){
08036 a[ncompstr*n+j]=dfds[j];
08037 }
08038
08039
08040
08041
08042
08043
08044 if (ncomphard>0){
08045
08046
08047
08048
08049 dgdsigmadq (ipp,im,sigt,q,dgdsdq);
08050
08051 mxm (d,dgdsdq,auxm2);
08052
08053 for (j=0;j<ncompstr;j++){
08054 for (k=0;k<ncomphard;k++){
08055 a[j*n+ncompstr+1+k]=auxm2[j][k]*dgamma;
08056 }
08057 }
08058
08059
08060
08061
08062
08063
08064 dhdsigma (ipp,im,sigt,q,dhds);
08065
08066
08067
08068
08069 for (j=0;j<ncomphard;j++){
08070 for (k=0;k<ncompstr;k++){
08071 a[(ncompstr+1+j)*n+k]=dhds[k][j]*dgamma;
08072 }
08073 }
08074
08075
08076
08077
08078 dfdqpar (ipp,im,sigt,q,dfdq);
08079
08080
08081
08082 for (j=0;j<ncomphard;j++){
08083 a[ncompstr*n+ncompstr+1+j]=dfdq[j];
08084 }
08085
08086
08087
08088
08089 dhdgamma (ipp,im,epsp,sigt,q,dhdc);
08090
08091
08092
08093 for (j=0;j<ncomphard;j++){
08094 a[(ncompstr+1+j)*n+ncompstr]=dhdc[j];
08095 }
08096
08097
08098
08099
08100 dhdqpar (ipp,im,ido,sigt,dgamma,q,dhdq);
08101
08102
08103
08104 for (j=0;j<ncomphard;j++){
08105 for (k=0;k<ncomphard;k++){
08106 a[(ncompstr+1+j)*n+ncompstr+1+k]=dhdq[j][k];
08107 }
08108 a[(ncompstr+1+j)*n+ncompstr+1+j]-=1.0;
08109 }
08110
08111 }
08112
08113
08114
08115
08116 for (j=0;j<ncompstr;j++){
08117 y[j]=sigtrial[j]-sig[j]-dgamma*auxv[j];
08118 }
08119
08120
08121 y[ncompstr]=0.0-f;
08122
08123
08124 for (j=0;j<ncomphard;j++){
08125
08126 }
08127
08128
08129
08130
08131
08132
08133
08134
08135
08136
08137
08138
08139
08140
08141
08142 gemp (a,x,y,n,1,1.0e-10,1);
08143
08144
08145
08146
08147
08148
08149
08150 norstr=0.0;
08151
08152 norstrincr=0.0;
08153
08154 for (j=0;j<ncompstr;j++){
08155 norstrincr+=x[j]*x[j];
08156 norstr+=sig[j]*sig[j];
08157 }
08158
08159 norconspar=x[ncompstr]*x[ncompstr];
08160
08161
08162 norhardpar=0.0;
08163
08164 norhardparincr=0.0;
08165
08166 for (j=0;j<ncomphard;j++){
08167 norhardpar+=q[j]*q[j];
08168 }
08169 for (j=ncompstr+1;j<ncompstr+1+ncomphard;j++){
08170 norhardparincr+=x[j]*x[j];
08171 }
08172
08173
08174
08175 for (j=0;j<ncompstr;j++){
08176 sig[j]+=x[j];
08177 auxv[j]=x[j];
08178 }
08179
08180
08181 dgamma+=x[ncompstr];
08182
08183
08184 for (j=0;j<ncomphard;j++){
08185 q[j]+=x[ncompstr+1+j];
08186 }
08187
08188
08189
08190
08191
08192
08193
08194
08195
08196
08197
08198
08199 mxv (c,auxv,auxv2);
08200
08201 for (j=0;j<ncompstr;j++){
08202 epsp[j]-=auxv2[j];
08203 }
08204
08205
08206
08207 vector_tensor (sig,sigt,ip[ipp].ssst,stress);
08208
08209
08210 f = yieldfunction (ipp,im,sigt,q);
08211
08212
08213 hardsoftfunction (ipp,im,sigt,q,hs);
08214
08215
08216
08217
08218
08219
08220
08221
08222
08223 j=0; k=0;
08224 if (norstr>zero){
08225 j++;
08226 if (norstrincr/norstr<err) k++;
08227 }
08228 if (norhardpar>zero){
08229 j++;
08230 if (norhardparincr/norhardpar<err) k++;
08231 }
08232 if (dgamma>zero){
08233 j++;
08234 if (norconspar/dgamma<err) k++;
08235 }
08236
08237
08238 if (j==k){
08239
08240
08241
08242
08243 break;
08244 }
08245 }
08246 }
08247
08248 fprintf (stdout,"\n\n parametr zpevneni %lf",q[0]);
08249 fprintf (stdout,"\n parametr konzistence %lf\n",gamma);
08250
08251
08252 gamma+=dgamma;
08253
08254
08255 storestress (0,ipp,sig);
08256
08257 delete [] a;
08258 delete [] x;
08259 delete [] y;
08260 }
08261
08262
08263
08264
08265
08266
08267
08268
08269
08270
08271
08272
08273
08274
08275 void mechmat::refresh_plast (long im, long ido)
08276 {
08277 if (plast == 0)
08278 for (long i=0;i<tnip;i++)
08279 if ( give_consparam(i,im,ido) ){
08280 plast = 1;
08281 return;
08282 }
08283 }
08284
08285
08286
08287
08288
08289
08290
08291
08292
08293
08294
08295
08296
08297
08298
08299
08300
08301
08302
08303
08304
08305
08306
08307
08308
08309
08310 void mechmat::damfuncpar(long ipp, long im, vector &eps, vector &kappa)
08311 {
08312 switch (ip[ipp].tm[im])
08313 {
08314 case scaldamage:
08315 {
08316 scdam[ip[ipp].idm[im]].damfuncpar(ipp, eps, kappa);
08317 break;
08318 }
08319 case scaldamagecc:
08320 {
08321 scdamcc[ip[ipp].idm[im]].damfuncpar(ipp, eps, kappa);
08322 break;
08323 }
08324 default:
08325 print_err("unknown scalar damage model is required",__FILE__,__LINE__,__func__);
08326 }
08327 return;
08328 }
08329
08330
08331
08332
08333
08334
08335
08336
08337
08338
08339
08340
08341
08342
08343
08344
08345
08346
08347
08348 double mechmat::damfunction(long ipp, long im, long ido,vector &kappa, vector &eps, vector &omegao)
08349 {
08350 double ret = 0.0;
08351
08352 switch (ip[ipp].tm[im])
08353 {
08354 case scaldamage:
08355 {
08356 ret = scdam[ip[ipp].idm[im]].damfunction(ipp,kappa,omegao);
08357 break;
08358 }
08359 case scaldamagecc:
08360 {
08361 ret = scdamcc[ip[ipp].idm[im]].damfunction(ipp,kappa,eps);
08362 break;
08363 }
08364 default:
08365 print_err("unknown scalar damage model is required",__FILE__,__LINE__,__func__);
08366 }
08367 return ret;
08368 }
08369
08370
08371
08372
08373
08374
08375
08376
08377
08378
08379
08380
08381
08382
08383 double mechmat::epsefunction(long ipp, long im)
08384 {
08385 double ret = 0.0;
08386
08387 switch (ip[ipp].tm[im])
08388 {
08389 case scaldamage:
08390 {
08391 ret = scdam[ip[ipp].idm[im]].epsefunction(ipp);
08392 break;
08393 }
08394 case scaldamagecc:
08395 {
08396 ret = scdamcc[ip[ipp].idm[im]].epsefunction();
08397 break;
08398 }
08399 default:
08400 print_err("unknown scalar damage model in function",__FILE__,__LINE__,__func__);
08401 }
08402 return ret;
08403 }
08404
08405
08406
08407
08408
08409
08410
08411
08412
08413
08414
08415
08416
08417
08418
08419
08420
08421
08422
08423
08424 double mechmat::scal_dam_sol (long ipp,long im,long ido,vector &eps,vector &kappa, vector &omegao, vector &sigma, matrix &d)
08425 {
08426 long nk = kappa.n;
08427
08428 double epse, omega;
08429 vector updkappa(nk);
08430
08431
08432
08433 if (nk > 1)
08434 {
08435 print_err("unsupported number of componetnts of vector kappa", __FILE__, __LINE__,__func__);
08436 return 0.0;
08437 }
08438
08439 if ((Mp->matmodel == local) || ((Mm->ip[ipp].hmt & 2) == 0))
08440
08441 {
08442 damfuncpar(ipp, im, eps, updkappa);
08443 if (updkappa[0] > kappa[0])
08444 {
08445 kappa[0] = updkappa[0];
08446 epse = epsefunction(ipp,im) ;
08447 if (kappa[0] > epse)
08448 omega = damfunction(ipp, im, ido, kappa, eps, omegao);
08449 else
08450 omega = 0.0;
08451 if (omega < omegao[0])
08452 omega = omegao[0];
08453 }
08454 else
08455 omega = omegao[0];
08456
08457 mxv(d, eps, sigma);
08458 cmulv(1.0-omega, sigma);
08459 return omega;
08460 }
08461
08462 if (Mp->matmodel == nonlocal)
08463
08464 {
08465 if (Mp->nonlocphase == 1)
08466
08467 {
08468 damfuncpar(ipp, im, eps, updkappa);
08469 kappa[0] = updkappa[0];
08470 fillv(0.0, sigma);
08471 return 0.0;
08472 }
08473 if (Mp->nonlocphase == 2)
08474
08475 {
08476 updkappa[0] = ip[ipp].nonloc[0];
08477 if (updkappa[0] > kappa[0])
08478 {
08479 kappa[0] = updkappa[0];
08480 epse = epsefunction(ipp, im);
08481 if (kappa[0] > epse)
08482 omega = damfunction(ipp, im, ido, kappa, eps, omegao);
08483 else
08484 omega = 0.0;
08485 if (omega < omegao[0])
08486 omega = omegao[0];
08487 }
08488 else
08489 omega = omegao[0];
08490
08491 mxv(d, eps, sigma);
08492 cmulv(1.0-omega, sigma);
08493 return omega;
08494 }
08495 }
08496
08497
08498 return 0.0;
08499 }
08500
08501
08502
08503
08504
08505
08506
08507
08508
08509
08510
08511
08512
08513
08514
08515
08516
08517
08518
08519
08520
08521
08522
08523
08524
08525
08526
08527
08528 void mechmat::givestressincr (long lcid,long ipp,long im,long ido,long fi,vector &sig)
08529 {
08530 long i;
08531
08532 if (im == 0);
08533 fillv(0.0, sig);
08534
08535 switch (ip[ipp].tm[im]){
08536 case effective_stress:{
08537 givestressincr (lcid, ipp, im+1, ido, fi, sig);
08538 break;
08539 }
08540 case shrinkagemat:{
08541 i=ip[ipp].idm[im];
08542 shmat[i].givestressincr (lcid, ipp, im, ido, fi, sig);
08543 break;
08544 }
08545 case creep_damage:{
08546 givestressincr (lcid, ipp, im+1, ido, fi, sig);
08547 break;
08548 }
08549
08550
08551
08552
08553
08554
08555 case viscoplasticity:{
08556 i=ip[ipp].idm[im];
08557 visplas[i].givestressincr (ipp,ido,fi,sig);
08558 break;
08559 }
08560 case creepb3:
08561 case creeprs:
08562 case creepdpl:{
08563 creep_givestressincr (ipp,im,ido,fi,sig);
08564 break;
08565 }
08566 case hypoplastmat:{
08567 i=ip[ipp].idm[im];
08568 hypopl[i].givestressincr(ipp, im, ido, sig);
08569 break;
08570 }
08571
08572
08573
08574
08575
08576
08577
08578
08579
08580
08581
08582
08583
08584
08585
08586
08587
08588
08589
08590
08591
08592
08593
08594
08595
08596
08597
08598
08599
08600
08601
08602
08603
08604 case elisomat:
08605 case elortomat:
08606 case simplas1d:
08607 case simvisc:
08608 case jflow:
08609 case jflow2:
08610 case microplaneM4:
08611 case microsimp:
08612 case microfibro:
08613 case mohrcoul:
08614 case mohrcoulparab:
08615 case boermaterial:
08616 case druckerprager:
08617 case doubledrprager:
08618 case druckerprager2:
08619 case chenplast:
08620 case modcamclaymat:
08621 case modcamclaycoupmat:
08622 case shefpl:
08623 case glasgowmechmat:
08624 case glasgowdamage:
08625 case scaldamage:
08626 case fixortodamage:
08627 case scaldamagecc:
08628 case ortodamage:
08629 case ortodamagerot:
08630 case anisodamage:
08631 case anisodamagerot:
08632 case graphm:
08633 case varelisomat:
08634 case geoelast:
08635 case creepeffym:
08636 case winklerpasternak:
08637 case nonlocdamgmat:
08638 case nonlocplastmat:
08639 case nonlocalmicroM4:
08640 case damage_plasticity:
08641 case cebfipcontmat:
08642 case contmat:
08643 break;
08644 default:{
08645 print_err("unknown material model is required",__FILE__,__LINE__,__func__);
08646 }
08647 }
08648 }
08649
08650
08651
08652
08653
08654
08655
08656
08657
08658
08659
08660
08661
08662
08663
08664
08665
08666
08667 void mechmat::stiff_deps_vispl (long ipp,long im,long ido,vector &sig,vector &q,double dt)
08668 {
08669 long i,ncomp,ncomph,ncompo;
08670
08671 ncomp=sig.n;
08672
08673 ncomph=q.n;
08674
08675
08676 ncompo = givencompother(ipp,im+1);
08677
08678 double f,g,dlambda,cs;
08679 vector epsp(ncomp),dq(ncomph),dh(ncomph),dfds(ncomp);
08680 matrix d(ncomp,ncomp),sigt(3,3),epspt(3,3), h(ncomph, ncomph);
08681
08682 vector_tensor (sig,sigt,ip[ipp].ssst,stress);
08683
08684 f = yieldfunction (ipp,im+2,sigt,q);
08685
08686 dfdsigma (ipp,im+2,sigt,q,dfds);
08687
08688 dfdqpar (ipp,im+2,sigt,q,dq);
08689
08690 switch (ip[ipp].tm[im+1]){
08691
08692
08693
08694
08695
08696 case simvisc:{
08697
08698 g = svis[ip[ipp].idm[im+1]].gfun (f);
08699 break;
08700 }
08701 case lemaitr:{
08702
08703 cs=ip[ipp].other[ido+3*ncomp];
08704
08705 g = lmtr[ip[ipp].idm[im+1]].gfun (f,cs);
08706 break;
08707 }
08708 default:{
08709 print_err("unknown material type is required",__FILE__,__LINE__,__func__);
08710 }
08711 }
08712
08713 dlambda = 0.0;
08714 if (g>0.0){
08715
08716 dlambda = dt*g;
08717
08718 for (i=0;i<ncomp;i++){
08719 epsp[i]=dlambda*dfds[i];
08720 }
08721
08722
08723 plasticmoduli (ipp, im+2, epsp, sigt, h);
08724 mxv (h,dq,dh);
08725 for (i=0;i<ncomph;i++){
08726 dh[i]=-dlambda*dh[i];
08727 }
08728 }
08729
08730
08731 switch (ip[ipp].tm[im+1]){
08732 case simvisc:{
08733
08734 for (i=0;i<ncomp;i++){
08735 ip[ipp].eqother[ido+ncompo+i]+=epsp[i];
08736 }
08737
08738 ip[ipp].eqother[ido+ncompo+ncomp]+=dlambda;
08739
08740 for (i=0;i<ncomph;i++){
08741 ip[ipp].eqother[ido+ncompo+ncomp+1+i]+=dh[i];
08742 }
08743
08744 break;
08745 }
08746
08747
08748
08749
08750
08751
08752
08753
08754
08755 case lemaitr:{
08756
08757 for (i=0;i<ncomp;i++){
08758 ip[ipp].eqother[ido+ncompo+i]+=epsp[i];
08759 }
08760
08761 vector_tensor (epsp,epspt,ip[ipp].ssst,strain);
08762 cumulstrain (epspt,cs);
08763 ip[ipp].eqother[ido+ncomp*2]=cs;
08764
08765
08766 ip[ipp].eqother[ido+ncompo+ncomp]+=dlambda;
08767
08768
08769 for (i=0;i<ncomph;i++){
08770 ip[ipp].eqother[ido+ncompo+ncomp+1+i]+=dh[i];
08771 }
08772 break;
08773 }
08774 default:{
08775 print_err("unknown material type is requiered",__FILE__,__LINE__,__func__);
08776 }
08777 }
08778
08779
08780 matstiff (d,ipp);
08781 mxv (d,epsp,sig);
08782
08783 }
08784
08785
08786
08787
08788
08789
08790
08791
08792
08793
08794
08795
08796
08797
08798 void mechmat::stiff_eps (long ipp)
08799 {
08800 long i,ncomp=ip[ipp].ncompstr;
08801 double time,
08802 complian,
08803 strength,
08804 shrink;
08805 vector eps(ncomp),sig(ncomp);
08806
08807 time=Mp->time/3600.0/24.0;
08808
08809 switch (ip[ipp].tm[0]){
08810 case aci:{
08811 aci78mod[ip[ipp].idm[0]].compliance (time,complian,strength,shrink);
08812 break;
08813 }
08814 case cebfip:{
08815 cebfip78mod[ip[ipp].idm[0]].compliance (time,complian,strength,shrink);
08816 break;
08817 }
08818 default:
08819 print_err("unknown material type is requiered",__FILE__,__LINE__,__func__);
08820 }
08821
08822
08823
08824 for (i=0;i<ncomp;i++){
08825 ip[ipp].strain[i]=ip[ipp].stress[i]*complian;
08826 }
08827 }
08828
08829
08830
08831
08832
08833
08834
08835
08836
08837
08838
08839
08840
08841
08842
08843
08844
08845
08846
08847
08848
08849
08850
08851 double mechmat::give_first_derivative (long ipp,double r)
08852 {
08853 double d;
08854
08855 switch (ip[ipp].tm[0]){
08856 case lenjonespot:{
08857 d = lenjon[ip[ipp].idm[0]].first_derivative (r);
08858 break;
08859 }
08860 default:{
08861 print_err("unknown material model is required",__FILE__,__LINE__,__func__);
08862 }
08863 }
08864
08865 return d;
08866 }
08867
08868
08869
08870
08871
08872
08873
08874
08875
08876
08877
08878
08879
08880 double mechmat::give_second_derivative (long ipp,double r)
08881 {
08882 double d;
08883
08884 switch (ip[ipp].tm[0]){
08885 case lenjonespot:{
08886 d = lenjon[ip[ipp].idm[0]].second_derivative (r);
08887 break;
08888 }
08889 default:{
08890 print_err("unknown material model is required",__FILE__,__LINE__,__func__);
08891 }
08892 }
08893
08894 return d;
08895 }
08896
08897
08898
08899
08900
08901
08902
08903
08904
08905
08906
08907
08908
08909
08910
08911
08912
08913
08914
08915
08916
08917
08918
08919
08920
08921
08922
08923
08924
08925
08926
08927
08928
08929
08930
08931
08932
08933
08934
08935
08936
08937
08938
08939
08940
08941
08942
08943
08944
08945
08946
08947
08948
08949
08950
08951
08952
08953
08954
08955
08956
08957
08958
08959
08960
08961
08962
08963
08964
08965
08966
08967
08968
08969
08970
08971
08972
08973
08974
08975
08976
08977
08978
08979
08980
08981
08982
08983
08984
08985
08986
08987
08988
08989
08990
08991
08992
08993
08994
08995
08996
08997
08998
08999
09000
09001
09002 void mechmat::save_intpoints_txt (FILE *aux, sel &selelems, sel *selother)
09003 {
09004 long i,j,k,n,nip,ir,ipp;
09005 sel selno;
09006 int prec = (int)Mp->hdbcont.prec;
09007
09008 selno.n=1;
09009 fprintf(aux, "\n");
09010
09011 if (nonmechq)
09012 {
09013 for (i=0; i<Mm->tnknmq; i++)
09014 fprintf(aux, "%ld ", Mm->nmqid[i]);
09015 fprintf(aux, "\n");
09016 }
09017
09018 for (i=0;i<Mt->ne;i++)
09019 {
09020 ipp = Mt->elements[i].ipp[0][0];
09021 nip = Mt->give_totnip(i);
09022 selelems.presence_id(i, ir);
09023 if (ir < 0)
09024 n = 0;
09025 else
09026 n = selother[ir].give_nselcomp(ip[ipp].ncompeqother);
09027 for (j=0; j<nip; j++)
09028 {
09029 fprintf (aux,"%ld %ld %ld %ld\n",ipp, Mb->nlc, ip[ipp].ncompstr, n);
09030 if (ir < 0)
09031 ip[ipp].save_data_txt(aux,Mb->nlc, selno);
09032 else
09033 ip[ipp].save_data_txt(aux,Mb->nlc, selother[ir]);
09034
09035
09036 if (nonmechq)
09037 {
09038 for (k=0; k<nnmq; k++)
09039 fprintf (aux,"%.*le ",prec, Mm->nonmechq[k][ipp]);
09040 fprintf(aux, "\n");
09041 }
09042 ipp++;
09043 }
09044 }
09045 }
09046
09047
09048
09049
09050
09051
09052
09053
09054
09055
09056
09057
09058
09059
09060
09061
09062
09063
09064
09065 void mechmat::save_intpoints_txt (long ni, sel &selelems, sel *selother)
09066 {
09067 long i,j,k,ir,n,nip,ipp;
09068 sel selno;
09069 selno.n=1;
09070 FILE *aux;
09071 char name[FNAMELEN+20];
09072 char emsg[FNAMELEN+100];
09073 int prec = (int)Mp->hdbcont.prec;
09074
09075 sprintf(name, "%s.%ld.strain.bac",Mp->hdbcont.hdbnames, ni);
09076 aux = fopen(name,"wt");
09077 if (aux==NULL)
09078 {
09079 sprintf(emsg, "cannot open backup file %s", name);
09080 print_err(emsg, __FILE__, __LINE__, __func__);
09081 abort();
09082 }
09083 for (i=0;i<tnip;i++)
09084 {
09085 fprintf (aux,"%ld %ld %ld\n",i, Mb->nlc, ip[i].ncompstr);
09086 for (j=0;j<Mb->nlc*ip[i].ncompstr;j++)
09087 fprintf (aux, "%.*le\n",prec , ip[i].strain[j]);
09088 }
09089 fclose(aux);
09090
09091 sprintf(name, "%s.%ld.stress.bac",Mp->hdbcont.hdbnames, ni);
09092 aux = fopen(name,"wt");
09093 if (aux==NULL)
09094 {
09095 sprintf(emsg, "cannot open backup file %s", name);
09096 print_err(emsg, __FILE__, __LINE__, __func__);
09097 abort();
09098 }
09099 for (i=0;i<tnip;i++)
09100 {
09101 fprintf (aux,"%ld %ld %ld\n",i, Mb->nlc, ip[i].ncompstr);
09102 for (j=0;j<Mb->nlc*ip[i].ncompstr;j++)
09103 fprintf (aux, "%.*le\n",prec , ip[i].stress[j]);
09104 }
09105 fclose(aux);
09106
09107 sprintf(name, "%s.%ld.other.bac",Mp->hdbcont.hdbnames, ni);
09108 aux = fopen(name,"wt");
09109 if (aux==NULL)
09110 {
09111 sprintf(emsg, "cannot open backup file %s", name);
09112 print_err(emsg, __FILE__, __LINE__, __func__);
09113 abort();
09114 }
09115 for (i=0;i<Mt->ne;i++)
09116 {
09117 ipp = Mt->elements[i].ipp[0][0];
09118 nip = Mt->give_totnip(i);
09119 selelems.presence_id(i, ir);
09120 if (ir < 0)
09121 n = 0;
09122 else
09123 n = selother[ir].give_nselcomp(ip[ipp].ncompeqother);
09124 for (j=0; j<nip; j++)
09125 {
09126 fprintf (aux,"%ld %ld\n",ipp, n);
09127 if (n==0)
09128 {
09129 ipp++;
09130 continue;
09131 }
09132 for (k=0;k<ip[ipp].ncompeqother;k++)
09133 {
09134 if (selother[ir].presence_id(k))
09135 fprintf (aux, "%.*le\n",prec , ip[ipp].eqother[k]);
09136 }
09137 ipp++;
09138 }
09139 }
09140 fclose(aux);
09141
09142
09143 if (nonmechq)
09144 {
09145 sprintf(name, "%s.%ld.trfquant.bac",Mp->hdbcont.hdbnames, ni);
09146 aux = fopen(name,"wt");
09147 if (aux==NULL)
09148 {
09149 sprintf(emsg, "cannot open backup file %s", name);
09150 print_err(emsg, __FILE__, __LINE__, __func__);
09151 abort();
09152 }
09153
09154 for (i=0; i<Mm->tnknmq; i++)
09155 fprintf(aux, "%ld ", Mm->nmqid[i]);
09156 fprintf(aux, "\n");
09157
09158 for (i=0;i<tnip;i++)
09159 {
09160 fprintf (aux,"%ld ",i);
09161 for (j=0; j<nnmq; j++)
09162 fprintf (aux,"%.*le ",prec , nonmechq[j][i]);
09163
09164 fprintf(aux, "\n");
09165 }
09166 fclose(aux);
09167 }
09168 }
09169
09170
09171
09172
09173
09174
09175
09176
09177
09178
09179
09180
09181
09182
09183
09184
09185
09186 void mechmat::save_intpoints_bin (FILE *aux, sel &selelems, sel *selother)
09187 {
09188 long i,j,k,n,nip,ir,ipp;
09189 sel selno;
09190
09191 selno.n=1;
09192
09193 if (nonmechq)
09194 fwrite (nmqid, sizeof(*nmqid), sizeof(nmqid), aux);
09195
09196 for (i=0;i<Mt->ne;i++)
09197 {
09198 ipp = Mt->elements[i].ipp[0][0];
09199 nip = Mt->give_totnip(i);
09200 selelems.presence_id(i, ir);
09201
09202 if (ir < 0)
09203 n = 0;
09204 else
09205 n = selother[ir].give_nselcomp(ip[ipp].ncompeqother);
09206 for (j=0; j<nip; j++)
09207 {
09208 fwrite(&ipp, sizeof(ipp), 1, aux);
09209 fwrite(&Mb->nlc, sizeof(Mb->nlc), 1, aux);
09210 fwrite(&ip[ipp].ncompstr, sizeof(ip[ipp].ncompstr), 1, aux);
09211 fwrite(&n, sizeof(n), 1, aux);
09212 if (ir < 0)
09213 ip[ipp].save_data_bin(aux,Mb->nlc, selno);
09214 else
09215 ip[ipp].save_data_bin(aux,Mb->nlc, selother[ir]);
09216
09217
09218 if (nonmechq)
09219 {
09220 for (k=0; k<nnmq; k++)
09221 fwrite (&Mm->nonmechq[k][ipp], sizeof(Mm->nonmechq[k][ipp]), 1, aux);
09222 }
09223 ipp++;
09224 }
09225 }
09226 }
09227
09228
09229
09230
09231
09232
09233
09234
09235
09236
09237
09238
09239
09240
09241
09242
09243
09244 void mechmat::save_intpoints_bin (long ni, sel &selelems, sel *selother)
09245 {
09246 long i,j,k,ir,n,nip,ipp;
09247 sel selno;
09248 selno.n=1;
09249 FILE *aux;
09250 char name[FNAMELEN+20];
09251 char emsg[FNAMELEN+100];
09252 double tmp;
09253
09254 sprintf(name, "%s.%ld.strain.bac",Mp->hdbcont.hdbnames, ni);
09255 aux = fopen(name,"wb");
09256 if (aux==NULL)
09257 {
09258 sprintf(emsg, "cannot open backup file %s", name);
09259 print_err(emsg, __FILE__, __LINE__, __func__);
09260 abort();
09261 }
09262 for (i=0;i<tnip;i++){
09263 fwrite(&i, sizeof(i), 1, aux);
09264 fwrite(&Mb->nlc, sizeof(Mb->nlc), 1, aux);
09265 fwrite(&ip[i].ncompstr, sizeof(ip[i].ncompstr), 1, aux);
09266 fwrite (ip[i].strain, sizeof(*ip[i].strain),Mb->nlc*ip[i].ncompstr, aux);
09267 }
09268 fclose(aux);
09269
09270 sprintf(name, "%s.%ld.stress.bac",Mp->hdbcont.hdbnames, ni);
09271 aux = fopen(name,"wb");
09272 if (aux==NULL)
09273 {
09274 sprintf(emsg, "cannot open backup file %s", name);
09275 print_err(emsg, __FILE__, __LINE__, __func__);
09276 abort();
09277 }
09278 for (i=0;i<tnip;i++)
09279 {
09280 fwrite(&i, sizeof(i), 1, aux);
09281 fwrite(&Mb->nlc, sizeof(Mb->nlc), 1, aux);
09282 fwrite(&ip[i].ncompstr, sizeof(ip[i].ncompstr), 1, aux);
09283 fwrite(ip[i].stress, sizeof(*ip[i].stress),Mb->nlc*ip[i].ncompstr, aux);
09284 }
09285 fclose(aux);
09286
09287 sprintf(name, "%s.%ld.other.bac",Mp->hdbcont.hdbnames, ni);
09288 aux = fopen(name,"wb");
09289 if (aux==NULL)
09290 {
09291 sprintf(emsg, "cannot open backup file %s", name);
09292 print_err(emsg, __FILE__, __LINE__, __func__);
09293 abort();
09294 }
09295 for (i=0;i<Mt->ne;i++)
09296 {
09297 ipp = Mt->elements[i].ipp[0][0];
09298 nip = Mt->give_totnip(i);
09299 selelems.presence_id(i, ir);
09300 if (ir < 0)
09301 n = 0;
09302 else
09303 n = selother[ir].give_nselcomp(ip[ipp].ncompeqother);
09304 for (j=0; j<nip; j++)
09305 {
09306 fwrite(&ipp, sizeof(ipp), 1, aux);
09307 fwrite(&n, sizeof(n), 1, aux);
09308 if (n==0)
09309 {
09310 ipp++;
09311 continue;
09312 }
09313 for (k=0;k<n;k++)
09314 {
09315 if (selother[ir].presence_id(k)){
09316 tmp = ip[ipp].eqother[k];
09317
09318 fwrite (&tmp, sizeof(tmp), 1, aux);
09319 }
09320 }
09321 ipp++;
09322 }
09323 }
09324 fclose(aux);
09325
09326 if (nonmechq)
09327 {
09328
09329 sprintf(name, "%s.%ld.trfquant.bac",Mp->hdbcont.hdbnames, ni);
09330 aux = fopen(name,"wb");
09331 if (aux==NULL)
09332 {
09333 sprintf(emsg, "cannot open backup file %s", name);
09334 print_err(emsg, __FILE__, __LINE__, __func__);
09335 abort();
09336 }
09337
09338 fwrite (nmqid, sizeof(*nmqid), sizeof(nmqid), aux);
09339
09340 for (ipp=0;ipp<tnip;ipp++)
09341 {
09342 fwrite(&ipp, sizeof(ipp), 1, aux);
09343 for (j=0; j<nnmq; j++)
09344 fwrite(&Mm->nonmechq[j][ipp], sizeof(Mm->nonmechq[j][ipp]), 1, aux);
09345 }
09346 fclose(aux);
09347 }
09348 }
09349
09350
09351
09352
09353
09354
09355
09356
09357
09358
09359
09360
09361
09362
09363
09364
09365
09366
09367
09368
09369 void mechmat::restore_intpoints_txt (FILE *aux, sel &selelems, sel *selother, long **selid)
09370 {
09371 long i,j,n,ir,ipp;
09372 long nlc, ncompstr;
09373 sel selno;
09374
09375 selno.n=1;
09376
09377 if (nonmechq)
09378 {
09379 for (i=0; i<Mm->tnknmq; i++)
09380 fscanf(aux, "%ld", &Mm->nmqid[i]);
09381 }
09382
09383 for (i=0;i<tnip;i++)
09384 {
09385 fscanf (aux,"%ld %ld %ld %ld",&ipp, &nlc, &ncompstr, &n);
09386 if ((ipp < 0) || (ipp >= tnip))
09387 {
09388 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09389 abort();
09390 }
09391 if (nlc != Mb->nlc)
09392 {
09393 print_err("invalid number of load cases", __FILE__, __LINE__, __func__);
09394 abort();
09395 }
09396 if (ncompstr != ip[ipp].ncompstr)
09397 {
09398 print_err("incompatible number of stress/strain components", __FILE__, __LINE__, __func__);
09399 abort();
09400 }
09401 selelems.presence_id(elip[ipp], ir);
09402 if (ir < 0)
09403 ip[ipp].restore_data_txt(aux, nlc, n, selno, NULL);
09404 else
09405 ip[ipp].restore_data_txt(aux, nlc, n, selother[ir], selid[ir]);
09406
09407
09408 if (nonmechq)
09409 {
09410 for (j=0; j<nnmq; j++)
09411 fscanf(aux,"%le", &Mm->nonmechq[j][ipp]);
09412 }
09413 }
09414 }
09415
09416
09417
09418
09419
09420
09421
09422
09423
09424
09425
09426
09427
09428
09429
09430
09431
09432
09433
09434
09435 void mechmat::restore_intpoints_txt (sel &selelems, sel *selother, long **selid)
09436 {
09437 long i,j,k,n,ir,ik,is,ipp;
09438 long nlc, ncompstr;
09439 double tmp;
09440 FILE *aux;
09441 char name[FNAMELEN+20];
09442 char emsg[FNAMELEN+100];
09443
09444
09445 sprintf(name, "%s.strain.bac",Mp->hdbcont.hdbnamer);
09446 aux = fopen(name,"rt");
09447 if (aux==NULL)
09448 {
09449 sprintf(emsg, "cannot open backup file %s", name);
09450 print_err(emsg, __FILE__, __LINE__, __func__);
09451 abort();
09452 }
09453 for (i=0;i<tnip;i++){
09454 fscanf (aux,"%ld %ld %ld",&ipp, &nlc, &ncompstr);
09455 if ((ipp < 0) || (ipp >= tnip))
09456 {
09457 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09458 abort();
09459 }
09460 if (nlc != Mb->nlc)
09461 {
09462 print_err("invalid number of load cases", __FILE__, __LINE__, __func__);
09463 abort();
09464 }
09465 if (ncompstr != ip[ipp].ncompstr)
09466 {
09467 print_err("incompatible number of stress/strain components", __FILE__, __LINE__, __func__);
09468 abort();
09469 }
09470 for (j=0; j<Mb->nlc*ip[i].ncompstr; j++)
09471 fscanf (aux, "%le", ip[i].strain+j);
09472 }
09473
09474 sprintf(name, "%s.stress.bac",Mp->hdbcont.hdbnamer);
09475 aux = fopen(name,"rt");
09476 if (aux==NULL)
09477 {
09478 sprintf(emsg, "cannot open backup file %s", name);
09479 print_err(emsg, __FILE__, __LINE__, __func__);
09480 abort();
09481 }
09482 for (i=0;i<tnip;i++){
09483 fscanf (aux,"%ld %ld %ld",&ipp, &nlc, &ncompstr);
09484 if ((ipp < 0) || (ipp >= tnip))
09485 {
09486 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09487 abort();
09488 }
09489 if (nlc != Mb->nlc)
09490 {
09491 print_err("invalid number of load cases", __FILE__, __LINE__, __func__);
09492 abort();
09493 }
09494 if (ncompstr != ip[ipp].ncompstr)
09495 {
09496 print_err("incompatible number of stress/strain components", __FILE__, __LINE__, __func__);
09497 abort();
09498 }
09499 for (j=0;j<nlc*ncompstr;j++)
09500 fscanf (aux, "%le", ip[i].stress+j);
09501 }
09502
09503 sprintf(name, "%s.other.bac",Mp->hdbcont.hdbnamer);
09504 aux = fopen(name,"rt");
09505 if (aux==NULL)
09506 {
09507 sprintf(emsg, "cannot open backup file %s", name);
09508 print_err(emsg, __FILE__, __LINE__, __func__);
09509 abort();
09510 }
09511 for (i=0;i<tnip;i++)
09512 {
09513 fscanf (aux,"%ld %ld",&ipp, &n);
09514 if ((ipp < 0) || (ipp >= tnip))
09515 {
09516 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09517 abort();
09518 }
09519 selelems.presence_id(elip[ipp], ir);
09520 for (k=0;k<n;k++)
09521 {
09522 fscanf (aux, "%le", &tmp);
09523 if (ir < 0)
09524 continue;
09525 if (selother[ir].presence_id(k,ik))
09526 {
09527 is = selid[ir][ik]+k-selother[ir].id1[ik];
09528 if (is >= ip[ipp].ncompeqother)
09529 print_err("invalid index for eqother restoring is required", __FILE__, __LINE__, __func__);
09530 else
09531 ip[ipp].eqother[is] = tmp;
09532 }
09533 }
09534 ipp++;
09535 }
09536 fclose(aux);
09537
09538
09539 if (nonmechq)
09540 {
09541 sprintf(name, "%s.trfquant.bac",Mp->hdbcont.hdbnamer);
09542 aux = fopen(name,"rt");
09543 if (aux==NULL)
09544 {
09545 sprintf(emsg, "cannot open backup file %s", name);
09546 print_err(emsg, __FILE__, __LINE__, __func__);
09547 abort();
09548 }
09549
09550
09551 for (i=0; i<Mm->tnknmq; i++)
09552 fscanf(aux, "%ld ", &Mm->nmqid[i]);
09553
09554 for (i=0;i<tnip;i++)
09555 {
09556 fscanf (aux, "%ld", &ipp);
09557 if ((ipp < 0) || (ipp >= tnip))
09558 {
09559 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09560 abort();
09561 }
09562 for (j=0; j<nnmq; j++)
09563 fscanf (aux,"%le", &nonmechq[j][ipp]);
09564 }
09565 fclose(aux);
09566 }
09567 }
09568
09569
09570
09571
09572
09573
09574
09575
09576
09577
09578
09579
09580
09581
09582
09583
09584
09585
09586
09587
09588 void mechmat::restore_intpoints_bin (FILE *aux, sel &selelems, sel *selother, long **selid)
09589 {
09590 long i, j, ipp, nlc, ncompstr, n;
09591 long ir;
09592 sel selno;
09593
09594 if (nonmechq)
09595 fread (nmqid, sizeof(*nmqid), sizeof(nmqid), aux);
09596
09597 for (i=0;i<tnip;i++){
09598 fread(&ipp, sizeof(ipp), 1, aux);
09599 fread(&nlc, sizeof(nlc), 1, aux);
09600 fread(&ncompstr, sizeof(ncompstr), 1, aux);
09601 fread(&n, sizeof(n), 1, aux);
09602 if ((ipp < 0) || (ipp >= tnip))
09603 {
09604 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09605 abort();
09606 }
09607 if (nlc != Mb->nlc)
09608 {
09609 print_err("invalid number of load cases", __FILE__, __LINE__, __func__);
09610 abort();
09611 }
09612 if (ncompstr != ip[ipp].ncompstr)
09613 {
09614 print_err("incompatible number of stress/strain components", __FILE__, __LINE__, __func__);
09615 abort();
09616 }
09617 selelems.presence_id(elip[ipp], ir);
09618 if (ir < 0)
09619 ip[ipp].restore_data_bin(aux, Mb->nlc, n, selno, NULL);
09620 else
09621 ip[ipp].restore_data_bin(aux, Mb->nlc, n, selother[ir], selid[ir]);
09622
09623
09624 if (nonmechq)
09625 {
09626 for (j=0; j<nnmq; j++)
09627 fread(&Mm->nonmechq[j][ipp], sizeof(Mm->nonmechq[j][ipp]), 1, aux);
09628 }
09629 }
09630 }
09631
09632
09633
09634
09635
09636
09637
09638
09639
09640
09641
09642
09643
09644
09645
09646
09647
09648
09649
09650 void mechmat::restore_intpoints_bin (sel &selelems, sel *selother, long **selid)
09651 {
09652 long i,k,n,ir,ik,is,ipp;
09653 long nlc, ncompstr;
09654 double tmp;
09655 FILE *aux;
09656 char name[FNAMELEN+20];
09657 char emsg[FNAMELEN+100];
09658
09659
09660 sprintf(name, "%s.strain.bac",Mp->hdbcont.hdbnamer);
09661 aux = fopen(name,"rb");
09662 if (aux==NULL)
09663 {
09664 sprintf(emsg, "cannot open backup file %s", name);
09665 print_err(emsg, __FILE__, __LINE__, __func__);
09666 abort();
09667 }
09668 for (i=0;i<tnip;i++){
09669 fread(&ipp, sizeof(ipp), 1, aux);
09670 fread(&nlc, sizeof(nlc), 1, aux);
09671 fread(&ncompstr, sizeof(ncompstr), 1, aux);
09672 if ((ipp < 0) || (ipp >= tnip))
09673 {
09674 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09675 abort();
09676 }
09677 if (nlc != Mb->nlc)
09678 {
09679 print_err("invalid number of load cases", __FILE__, __LINE__, __func__);
09680 abort();
09681 }
09682 if (ncompstr != ip[ipp].ncompstr)
09683 {
09684 print_err("incompatible number of stress/strain components", __FILE__, __LINE__, __func__);
09685 abort();
09686 }
09687 fread (ip[ipp].strain, sizeof(*ip[ipp].strain),Mb->nlc*ip[ipp].ncompstr, aux);
09688 }
09689
09690 sprintf(name, "%s.stress.bac",Mp->hdbcont.hdbnamer);
09691 aux = fopen(name,"rb");
09692 if (aux==NULL)
09693 {
09694 sprintf(emsg, "cannot open backup file %s", name);
09695 print_err(emsg, __FILE__, __LINE__, __func__);
09696 abort();
09697 }
09698 for (i=0;i<tnip;i++){
09699 fread(&ipp, sizeof(ipp), 1, aux);
09700 fread(&nlc, sizeof(nlc), 1, aux);
09701 fread(&ncompstr, sizeof(ncompstr), 1, aux);
09702 if ((ipp < 0) || (ipp >= tnip))
09703 {
09704 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09705 abort();
09706 }
09707 if (nlc != Mb->nlc)
09708 {
09709 print_err("invalid number of load cases", __FILE__, __LINE__, __func__);
09710 abort();
09711 }
09712 if (ncompstr != ip[ipp].ncompstr)
09713 {
09714 print_err("incompatible number of stress/strain components", __FILE__, __LINE__, __func__);
09715 abort();
09716 }
09717 fread (ip[ipp].stress, sizeof(*ip[ipp].stress),Mb->nlc*ip[ipp].ncompstr, aux);
09718 }
09719
09720 sprintf(name, "%s.other.bac",Mp->hdbcont.hdbnamer);
09721 aux = fopen(name,"rb");
09722 if (aux==NULL)
09723 {
09724 sprintf(emsg, "cannot open backup file %s", name);
09725 print_err(emsg, __FILE__, __LINE__, __func__);
09726 abort();
09727 }
09728 for (i=0;i<tnip;i++){
09729 fread (&ipp, sizeof(ipp), 1, aux);
09730 fread (&n, sizeof(n), 1, aux);
09731 if ((ipp < 0) || (ipp >= tnip))
09732 {
09733 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09734 abort();
09735 }
09736 selelems.presence_id(elip[ipp], ir);
09737 for (k=0;k<n;k++)
09738 {
09739 fread (&tmp, sizeof(tmp), 1, aux);
09740 if (ir < 0)
09741 continue;
09742 if (selother[ir].presence_id(k,ik))
09743 {
09744 is = selid[ir][ik]+k-selother[ir].id1[ik];
09745 if (is >= ip[ipp].ncompeqother)
09746 print_err("invalid index for eqother restoring is required", __FILE__, __LINE__, __func__);
09747 else
09748 ip[ipp].eqother[is] = tmp;
09749 }
09750 }
09751 }
09752 fclose(aux);
09753
09754 if (nonmechq)
09755 {
09756 sprintf(name, "%s.trfquant.bac",Mp->hdbcont.hdbnamer);
09757 aux = fopen(name,"rb");
09758 if (aux==NULL)
09759 {
09760 sprintf(emsg, "cannot open backup file %s", name);
09761 print_err(emsg, __FILE__, __LINE__, __func__);
09762 abort();
09763 }
09764
09765
09766 fread (nmqid, sizeof(*nmqid), sizeof(nmqid), aux);
09767
09768 for (i=0;i<tnip;i++)
09769 {
09770 fread (&ipp, sizeof(ipp), 1, aux);
09771 if ((ipp < 0) || (ipp >= tnip))
09772 {
09773 print_err("invalid integration point number", __FILE__, __LINE__, __func__);
09774 abort();
09775 }
09776 for (k=0; k<nnmq; k++)
09777 fread(&Mm->nonmechq[k][ipp], sizeof(Mm->nonmechq[k][ipp]), 1, aux);
09778 }
09779 fclose(aux);
09780 }
09781 }