00001 #include "fixortodam.h"
00002 #include "global.h"
00003 #include "sequent.h"
00004 #include "intpoints.h"
00005 #include "matrix.h"
00006 #include "vecttens.h"
00007
00008
00009 fixortodam::fixortodam (void)
00010 {
00011 x1t = x2t = x3t = 0.0;
00012 allocv(3, xeqtf);
00013 allocv(3, xeqt0);
00014 }
00015
00016
00017
00018 fixortodam::~fixortodam (void)
00019 {
00020 }
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 void fixortodam::read (XFILE *in)
00035 {
00036 xfscanf(in, "%k %le", "x1t", &x1t);
00037 xfscanf(in, "%k %le", "x2t", &x2t);
00038 xfscanf(in, "%k %le", "x3t", &x3t);
00039
00040 xfscanf(in, "%k %le", "xeq1tf", &xeqtf(0));
00041 xfscanf(in, "%k %le", "xeq2tf", &xeqtf(1));
00042 xfscanf(in, "%k %le", "xeq3tf", &xeqtf(2));
00043 }
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void fixortodam::print (FILE *out)
00058 {
00059 fprintf(out, "%le ", x1t);
00060 fprintf(out, "%le ", x2t);
00061 fprintf(out, "%le ", x3t);
00062
00063 fprintf(out, "%le ", xeqtf(0));
00064 fprintf(out, "%le ", xeqtf(1));
00065 fprintf(out, "%le ", xeqtf(2));
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078 void fixortodam::initval(long ipp, long ido)
00079 {
00080 matrix d(6,6);
00081
00082
00083
00084 Mm->loc_elmatstiff (d,ipp);
00085
00086
00087
00088
00089
00090 for(long i=0; i<3; i++)
00091 Mm->ip[ipp].eqother[ido+3+i] = 2*xeqtf(i);
00092
00093
00094
00095
00096 }
00097
00098
00099
00100 void fixortodam::matstiff (matrix &d,long ipp,long ido)
00101 {
00102 matrix dl(6,6), dg(6,6);
00103 matrix tmat(3,3);
00104
00105 switch (Mp->nlman->stmat)
00106 {
00107 case initial_stiff:
00108
00109
00110 Mm->loc_elmatstiff (dl,ipp);
00111 break;
00112 case tangent_stiff:
00113 tmatstiff(dl, ipp, ido);
00114 break;
00115 default:
00116 print_err("unknown type of stifness matrix is required.", __FILE__, __LINE__, __func__);
00117 }
00118
00119
00120 Mm->loc_transf_mat(tmat, ipp);
00121
00122
00123
00124 lg_tens4transf(dg, dl, tmat, Mm->ip[ipp].ssst);
00125
00126 tensor4_matrix(d, dg, Mm->ip[ipp].ssst);
00127 }
00128
00129
00130
00131 void fixortodam::tmatstiff (matrix &d,long ipp,long ido)
00132 {
00133 long i;
00134 vector omega(3);
00135
00136
00137 for (i=0; i<3; i++)
00138 omega(i) = Mm->ip[ipp].eqother[ido+i];
00139
00140 secstiffmat(d, ipp, omega);
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 void fixortodam::compute_eqdispl(long ipp, matrix &epst, vector &xeq)
00154 {
00155 long i;
00156 long n = Mm->ip[ipp].ncompstr;
00157 double l;
00158
00159
00160
00161
00162
00163
00164
00165 long eid = Mm->elip[ipp];
00166 switch (Mt->give_dimension(eid))
00167 {
00168 case 1:
00169 l = Mt->give_length(eid);
00170 break;
00171 case 2:
00172 l = sqrt(Mt->give_area(eid));
00173 break;
00174 case 3:
00175 l = pow(Mt->give_volume(eid), 1.0/3.0);
00176 break;
00177 default:
00178 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00179 }
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 for (i=0; i<3; i++)
00197 xeq(i)=l*fabs(epst(i,i));
00198
00199 }
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 void fixortodam::compute_dam (vector &xeq, vector &xeqt0, vector &omegao, vector &omega)
00214 {
00215 long i;
00216
00217 for (i=0; i<3; i++)
00218 {
00219 if (xeq(i) >= xeqtf(i))
00220 {
00221 omega(i) = 1.0;
00222 continue;
00223 }
00224
00225 if (xeq(i) <= xeqt0(i))
00226 {
00227 omega(i) = 0.0;
00228 continue;
00229 }
00230 omega(i)=xeqtf(i)*(xeq(i)-xeqt0(i))/(xeq(i)*(xeqtf(i)-xeqt0(i)));
00231 if (omega(i)<omegao(i))
00232 omega(i)=omegao(i);
00233 if (omega(i)<0)
00234 omega(i)=0;
00235 }
00236 }
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 void fixortodam::secstiffmat(matrix &d, long ipp, vector &omega)
00251 {
00252
00253
00254 Mm->loc_elmatstiff (d,ipp);
00255
00256 d(0,0)=pow(1-omega(0),2)*d(0,0);
00257 d(0,1)=(1-omega(0))*(1-omega(1))*d(0,1);
00258 d(0,2)=(1-omega(0))*(1-omega(2))*d(0,2);
00259 d(1,0)=d(0,1);
00260 d(1,1)=pow(1-omega(1),2)*d(1,1);
00261 d(1,2)=(1-omega(1))*(1-omega(2))*d(1,2);
00262 d(2,0)=d(0,2);
00263 d(2,1)=d(1,2);
00264 d(2,2)=pow(1-omega(2),2)*d(2,2);
00265 if (omega(0) == 1.0 && omega(1) == 1.0)
00266 d(3, 3) = 0;
00267 else
00268 d(3, 3) = pow((2 * (1 - omega(0))*(1 - omega(1))) / (2 - omega(0) - omega(1)), 2)*d(3,3);
00269 if (omega(0) == 1.0 && omega(2) == 1.0)
00270 d(4, 4) = 0;
00271 else
00272 d(4, 4) = pow((2 * (1 - omega(0))*(1 - omega(2))) / (2 - omega(0) - omega(2)), 2)*d(4,4);
00273 if (omega(1) == 1.0 && omega(2) == 1.0)
00274 d(5, 5) = 0;
00275 else
00276 d(5, 5) = pow((2 * (1 - omega(1))*(1 - omega(2))) / (2 - omega(1) - omega(2)), 2)*d(5,5);
00277
00278
00279
00280 }
00281
00282
00283
00284 void fixortodam::nlstresses (long ipp, long im, long ido)
00285 {
00286 long i;
00287 long ncompstr=Mm->ip[ipp].ncompstr;
00288
00289 vector omega(3), omegao(3), xeq(3);
00290 vector epsg(ncompstr), eps(ncompstr), sig(ncompstr), sig_g(ncompstr);
00291 matrix epst(3,3), sigt(3,3), d(6,6);
00292 matrix tmat(3,3);
00293
00294
00295
00296 vector xeqt0;
00297
00298
00299 for (i=0; i<ncompstr; i++)
00300 epsg(i) = Mm->ip[ipp].strain[i];
00301
00302
00303 Mm->loc_transf_mat(tmat, ipp);
00304
00305
00306 vector_tensor (epsg,epst,Mm->ip[ipp].ssst,strain);
00307 glmatrixtransf(epst, tmat);
00308 tensor_vector(eps,epst,Mm->ip[ipp].ssst,strain);
00309
00310
00311
00312
00313
00314 for (i=0; i<3; i++)
00315 {
00316 omegao(i) = Mm->ip[ipp].eqother[ido+i];
00317 xeqt0(i) = Mm->ip[ipp].eqother[ido+3+i];
00318 }
00319
00320
00321 compute_eqdispl(ipp, epst, xeq);
00322
00323
00324 compute_dam(xeq, xeqt0, omegao, omega);
00325
00326
00327 secstiffmat(d, ipp, omega);
00328
00329
00330
00331 mxv(d, eps, sig);
00332
00333 if (sig[0]>x1t)
00334 xeqt0[0] = eps[0];
00335 if (sig[1]>x2t)
00336 xeqt0[1] = eps[1];
00337 if (sig[2]>x3t)
00338 xeqt0[2] = eps[2];
00339
00340
00341 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00342 lgmatrixtransf(sigt, tmat);
00343 tensor_vector(sig_g,sigt,Mm->ip[ipp].ssst,stress);
00344
00345
00346 for (i=0; i<ncompstr; i++)
00347 Mm->ip[ipp].stress[i] = sig_g[i];
00348
00349
00350
00351 for (i=0; i<3; i++)
00352 Mm->ip[ipp].other[ido+i] = omega(i);
00353 for (i=0; i<3; i++)
00354 Mm->ip[ipp].other[ido+3+i] = xeqt0(i);
00355 }
00356
00357
00358
00359 void fixortodam::updateval (long ipp,long im,long ido)
00360 {
00361 long i;
00362
00363 for (i=0; i<6; i++)
00364 Mm->ip[ipp].eqother[ido+i] = Mm->ip[ipp].other[ido+i];
00365 }