00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include "hypoplast.h"
00005 #include "global.h"
00006 #include "intpoints.h"
00007 #include "vecttens.h"
00008
00009 #include "iotools.h"
00010 #include "matrix.h"
00011 #include "vector.h"
00012
00013 #ifdef FCC_DEF
00014 extern "C" void umatunsat_(double &, double& , double &, double &,
00015 double &, double &, double &, double &,
00016 double &, double &, double &, double &,
00017 double &, double &, double &, double &,
00018 double &, double &, char &, int &,
00019 int &, int &, int &, double &,
00020 int &, double &, double &, double &,
00021 double &, double &, double &, int &,
00022 int &, int &, int &, int &,
00023 int &, double &);
00024 #else
00025 void umatunsat_(double &, double& , double &, double &,
00026 double &, double &, double &, double &,
00027 double &, double &, double &, double &,
00028 double &, double &, double &, double &,
00029 double &, double &, char &, int &,
00030 int &, int &, int &, double &,
00031 int &, double &, double &, double &,
00032 double &, double &, double &, int &,
00033 int &, int &, int &, int &,
00034 int &, double &)
00035 {
00036 print_err("Hypoplastic material from ABAQUS was not compiled\n"
00037 "(Fortran compiler was not found)", __FILE__, __LINE__, __func__);
00038 }
00039 #endif
00040
00041
00042 hypoplast::hypoplast()
00043 {
00044 phi = p_t = lam_star = kap_star = n_star = 0.0;
00045 rr = n = l = m = s_e0 = e_0;
00046 nprops = 11;
00047 nstatv = 6;
00048 pr_suc_fl = no;
00049 props = new double[nprops];
00050 memset(props, 0, sizeof(*props)*nprops);
00051 }
00052
00053
00054
00055 hypoplast::~hypoplast()
00056 {
00057 }
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 void hypoplast::read(XFILE *in)
00072 {
00073 xfscanf(in, "%k %le %k %le %k %le %k %le %k %le",
00074 "phi", &phi, "p_t", &p_t, "lam_star", &lam_star,
00075 "kap_star", &kap_star, "N_star", &n_star);
00076 xfscanf(in, "%k %le %k %le %k %le %k %le %k %le %k %le",
00077 "rr", &rr, "n", &n, "l", &l, "m", &m,
00078 "s_e0", &s_e0, "e_0", &e_0);
00079 xfscanf(in, "%k%m", "presc_suction", &answertype_kwdset, &pr_suc_fl);
00080 if (pr_suc_fl == yes)
00081 suc_fn.read(in);
00082
00083 props[0] = phi;
00084 props[1] = p_t;
00085 props[2] = lam_star;
00086 props[3] = kap_star;
00087 props[4] = n_star;
00088 props[5] = rr;
00089 props[6] = n;
00090 props[7] = l;
00091 props[8] = m;
00092 props[9] = s_e0 ;
00093 props[10] = e_0;
00094 }
00095
00096
00097
00098 void hypoplast::initval (long ipp,long ido)
00099 {
00100 long ncomp = Mm->ip[ipp].ncompstr;
00101 double *statev = Mm->ip[ipp].eqother+ido+2*ncomp;
00102
00103 if (statev[0] == 0.0)
00104 {
00105 statev[0] = Mm->ic[ipp][0];
00106
00107
00108 if ((Mp->eigstrains == 4) || (Mp->eigstrains == 5))
00109 {
00110 for (long i=0; i<ncomp; i++)
00111 Mm->ip[ipp].eqother[ido+ncomp+i] = -Mm->eigstresses[ipp][i];
00112 }
00113 }
00114
00115
00116
00117
00118
00119 if (pr_suc_fl == yes)
00120 statev[1] = suc_fn.getval(Mp->time);
00121 else
00122 statev[1] = 200.0;
00123
00124 statev[2] = 0.0;
00125 }
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 void hypoplast::matstiff (matrix &d,long ipp,long ido)
00141 {
00142 long i;
00143 int noel = int(Mm->elip[ipp]);
00144 int npt = int(ipp);
00145 int ndi = 3;
00146 int nsh = 3;
00147 int ntens = ndi+nsh;
00148
00149
00150 double ddum=0.0;
00151 double ddum3[3]={0.0, 0.0, 0.0};
00152 double ddum6[6]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
00153 double dfgrd[3*3]={1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
00154 int idum;
00155 char cmname []="unsatm";
00156
00157
00158 double unsatvar[5]={0.0, 0.0, 0.0, 0.0, 0.0};
00159
00160
00161 long ncomp = Mm->ip[ipp].ncompstr;
00162 strastrestate ssst = Mm->ip[ipp].ssst;
00163 double *statev = Mm->ip[ipp].other+ido+2*ncomp;
00164 double aux, time, dtime;
00165 vector sig(6), eps(6);
00166 matrix dd(6,6);
00167
00168
00169 for (i=ido+2*ncomp; i<ido+2*ncomp+nstatv; i++)
00170 Mm->ip[ipp].other[i] = Mm->ip[ipp].eqother[i];
00171
00172
00173
00174
00175
00176
00177 if (pr_suc_fl == yes)
00178 {
00179 unsatvar[0] = Mm->ip[ipp].eqother[ido+2*ncomp+1];
00180 unsatvar[1] = suc_fn.getval(Mp->time)-unsatvar[0];
00181 }
00182 else
00183 {
00184 unsatvar[0] = 200.0;
00185 unsatvar[1] = 0.0;
00186 }
00187 unsatvar[2]= 0.0;
00188
00189 statev[1]=unsatvar[0];
00190 statev[2]=unsatvar[2];
00191
00192
00193 for (i=0; i<ncomp; i++)
00194 Mm->ip[ipp].stress[i] = Mm->ip[ipp].eqother[ido+ncomp+i];
00195
00196
00197 give_full_vector(sig.a, Mm->ip[ipp].stress, ssst);
00198 give_full_vector(eps.a, Mm->ip[ipp].strain, ssst);
00199
00200
00201
00202
00203
00204
00205 aux = sig[3];
00206 sig[3] = sig[5];
00207 sig[5] = aux;
00208
00209
00210
00211
00212 aux = eps[3];
00213 eps[3] = eps[5];
00214 eps[5] = aux;
00215
00216 if (Mp->timecon.tct == notct)
00217 {
00218 dtime = Mp->dlambda;
00219 time = Mp->lambda;
00220 }
00221 else
00222 {
00223 dtime = Mp->dtime;
00224 time = Mp->time;
00225 }
00226
00227 umatunsat_(*sig.a, *statev, *dd.a, ddum, ddum,
00228 ddum, ddum, *ddum6, *ddum6, ddum,
00229 *eps.a, *ddum6, time, dtime, ddum,
00230 ddum, *ddum6, *ddum6, *cmname, ndi,
00231 nsh, ntens, nstatv, *props, nprops,
00232 *ddum3, *dfgrd, ddum, ddum, *dfgrd,
00233 *dfgrd, noel, npt, idum, idum,
00234 (int&)Mp->istep, (int&)Mp->jstep, *unsatvar);
00235
00236
00237
00238
00239 for (i=0; i<5; i++)
00240 {
00241 if (i == 3)
00242 continue;
00243 aux = dd[i][3];
00244 dd[i][3] = dd[i][5];
00245 dd[i][5] = aux;
00246 }
00247
00248 for (i=0; i<5; i++)
00249 {
00250 if (i == 3)
00251 continue;
00252 aux = dd[3][i];
00253 dd[3][i] = dd[5][i];
00254 dd[5][i] = aux;
00255 }
00256
00257 aux = dd[3][3];
00258 dd[3][3] = dd[5][5];
00259 dd[5][5] = aux;
00260 aux = dd[3][5];
00261 dd[3][5] = dd[5][3];
00262 dd[5][3] = aux;
00263
00264
00265 tensor4_matrix (d, dd, ssst);
00266 }
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283 void hypoplast::givestressincr (long ipp, long im, long ido, vector &sig)
00284 {
00285 long i;
00286
00287
00288 long ncomp = Mm->ip[ipp].ncompstr;
00289 double ds;
00290 matrix sigt(3,3);
00291
00292 if (pr_suc_fl == yes)
00293 ds = suc_fn.getval(Mp->time)-Mm->ip[ipp].eqother[ido+2*ncomp+1];
00294 else
00295 ds = 0.0;
00296
00297
00298 for (i=0;i<3;i++)
00299 sigt[i][i] += ds;
00300
00301
00302 tensor_vector(sig, sigt, Mm->ip[ipp].ssst,stress);
00303 }
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319 void hypoplast::nlstresses (long ipp, long im, long ido)
00320 {
00321 long i;
00322 int noel = int(Mm->elip[ipp]);
00323 int npt = int(ipp);
00324 int ndi = 3;
00325 int nsh = 3;
00326 int ntens = ndi+nsh;
00327
00328
00329 double ddum=0.0;
00330 double ddum3[3]={0.0, 0.0, 0.0};
00331 double ddum6[6]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
00332 double ddum36[6*6];
00333 memset(ddum36, 0, sizeof(*ddum36)*36);
00334 double dfgrd[3*3]={1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
00335 int idum;
00336 char cmname []="unsatm";
00337
00338
00339 double unsatvar[5]={0.0, 0.0, 0.0, 0.0, 0.0};
00340
00341
00342 long ncomp = Mm->ip[ipp].ncompstr;
00343 strastrestate ssst = Mm->ip[ipp].ssst;
00344 double *statev = Mm->ip[ipp].other+ido+2*ncomp;
00345 double aux, time, dtime;
00346 double pnewdt = 1.0;
00347 vector sig(6), eps(6), deps(6);
00348 matrix dd(6,6);
00349
00350
00351
00352
00353 for (i=ido+2*ncomp; i<ido+2*ncomp+nstatv; i++)
00354 Mm->ip[ipp].other[i] = Mm->ip[ipp].eqother[i];
00355
00356
00357
00358
00359
00360
00361
00362 if (pr_suc_fl == yes)
00363 {
00364 unsatvar[0] = Mm->ip[ipp].eqother[ido+2*ncomp+1];
00365 unsatvar[1] = suc_fn.getval(Mp->time)-unsatvar[0];
00366 }
00367 else
00368 {
00369 unsatvar[0] = 200.0;
00370 unsatvar[1] = 0.0;
00371 }
00372 unsatvar[2]= 0.0;
00373
00374 statev[1]=unsatvar[0];
00375 statev[2]=unsatvar[2];
00376
00377
00378
00379 for (i=0; i<ncomp; i++)
00380 eps[i] = Mm->ip[ipp].strain[i] - Mm->ip[ipp].eqother[ido+i];
00381
00382
00383 for (i=0; i<ncomp; i++)
00384 Mm->ip[ipp].stress[i] = Mm->ip[ipp].eqother[ido+ncomp+i];
00385
00386
00387 give_full_vector(sig.a, Mm->ip[ipp].stress, ssst);
00388 give_full_vector(deps.a, eps.a, ssst);
00389 give_full_vector(eps.a, Mm->ip[ipp].strain, ssst);
00390
00391
00392
00393
00394
00395 aux = sig[3];
00396 sig[3] = sig[5];
00397 sig[5] = aux;
00398
00399
00400
00401
00402 aux = eps[3];
00403 eps[3] = eps[5];
00404 eps[5] = aux;
00405
00406
00407
00408
00409 aux = deps[3];
00410 deps[3] = deps[5];
00411 deps[5] = aux;
00412
00413 if (Mp->timecon.tct == notct)
00414 {
00415 dtime = Mp->dlambda;
00416 time = Mp->lambda;
00417 }
00418 else
00419 {
00420 dtime = Mp->dtime;
00421 time = Mp->time;
00422 }
00423
00424 umatunsat_(*sig.a, *statev, *ddum36, ddum, ddum,
00425 ddum, ddum, *ddum6, *ddum6, ddum,
00426 *eps.a, *deps.a, time, dtime, ddum,
00427 ddum, *ddum6, *ddum6, *cmname, ndi,
00428 nsh, ntens, nstatv, *props, nprops,
00429 *ddum3, *dfgrd, pnewdt, ddum, *dfgrd,
00430 *dfgrd, noel, npt, idum, idum,
00431 (int&)Mp->istep, (int&)Mp->jstep, *unsatvar);
00432
00433
00434
00435
00436 aux = sig[3];
00437 sig[3] = sig[5];
00438 sig[5] = aux;
00439
00440
00441 give_red_vector(sig.a, Mm->ip[ipp].stress, ssst);
00442
00443
00444 Mm->storeother(ipp, ido, ncomp, Mm->ip[ipp].strain);
00445
00446 Mm->storeother(ipp, ido+ncomp, ncomp, Mm->ip[ipp].stress);
00447
00448 Mm->ip[ipp].other[ido+2*ncomp+1] = unsatvar[0];
00449
00450 Mm->ip[ipp].other[ido+2*ncomp+nstatv]= pnewdt;
00451 }
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467 double hypoplast::dstep_red(long ipp, long im, long ido)
00468 {
00469 long ncomp = Mm->ip[ipp].ncompstr;
00470 return (Mm->ip[ipp].other[ido+2*ncomp+nstatv]);
00471
00472 }
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488 void hypoplast::updateval (long ipp,long im,long ido)
00489 {
00490 long i;
00491 long n = Mm->givencompeqother(ipp,im);
00492
00493 for (i=0;i<n;i++)
00494 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00495 }
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509 void hypoplast::give_reqnmq(long *anmq)
00510 {
00511 anmq[saturation_degree-1] = 1;
00512 anmq[suction-1] = 1;
00513 }
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528 double hypoplast::give_virgporosity(long ipp, long ido)
00529 {
00530
00531 return exp(n_star)-1.0;
00532 }
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547 double hypoplast::give_iniporosity(long ipp, long ido)
00548 {
00549 return Mm->ic[ipp][0];
00550 }
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564 void hypoplast::changeparam (atsel &atm,vector &val)
00565 {
00566 long i;
00567
00568 for (i=0;i<atm.num;i++){
00569 switch (atm.atrib[i]){
00570 case 0:{
00571 phi=val[i];
00572 break;
00573 }
00574 case 1:{
00575 p_t=val[i];
00576 break;
00577 }
00578 case 2:{
00579 lam_star=val[i];
00580 break;
00581 }
00582 case 3:{
00583 kap_star=val[i];
00584 break;
00585 }
00586 case 4:{
00587 n_star=val[i];
00588 break;
00589 }
00590 case 5:{
00591 rr=val[i];
00592 break;
00593 }
00594 case 6:{
00595 n=val[i];
00596 break;
00597 }
00598 case 7:{
00599 l=val[i];
00600 break;
00601 }
00602 case 8:{
00603 m=val[i];
00604 break;
00605 }
00606 case 9:{
00607 s_e0=val[i];
00608 break;
00609 }
00610 case 10:{
00611 e_0=val[i];
00612 break;
00613 }
00614 default:{
00615 fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
00616 }
00617 }
00618 }
00619 }