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