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 xeqt0(0) = x1t/d(0,0);
00087 xeqt0(1) = x2t/d(1,1);
00088 xeqt0(2) = x3t/d(2,2);
00089 }
00090
00091
00092
00093 void fixortodam::matstiff (matrix &d,long ipp,long ido)
00094 {
00095 matrix dl(6,6), dg(6,6);
00096 matrix tmat(3,3);
00097
00098 switch (Mp->nlman->stmat)
00099 {
00100 case initial_stiff:
00101
00102
00103 Mm->loc_elmatstiff (dl,ipp);
00104 break;
00105 case tangent_stiff:
00106 tmatstiff(dl, ipp, ido);
00107 break;
00108 default:
00109 print_err("unknown type of stifness matrix is required.", __FILE__, __LINE__, __func__);
00110 }
00111
00112
00113 Mm->loc_transf_mat(tmat, ipp);
00114
00115
00116
00117 lg_tens4transf(dg, dl, tmat, Mm->ip[ipp].ssst);
00118
00119 tensor4_matrix(d, dg, Mm->ip[ipp].ssst);
00120 }
00121
00122
00123
00124 void fixortodam::tmatstiff (matrix &d,long ipp,long ido)
00125 {
00126 long i;
00127 vector omega(3);
00128
00129
00130 for (i=0; i<3; i++)
00131 omega(i) = Mm->ip[ipp].eqother[ido+i];
00132
00133 secstiffmat(d, ipp, omega);
00134 }
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 void fixortodam::compute_eqdispl(long ipp, matrix &epst, vector &xeq)
00147 {
00148 long i;
00149 long n = Mm->ip[ipp].ncompstr;
00150 double l;
00151 matrix dl(6,6);
00152 vector eps(n);
00153 vector sig(n);
00154 matrix sigt(3,3);
00155
00156
00157
00158 Mm->loc_elmatstiff (dl,ipp);
00159
00160 long eid = Mm->elip[ipp];
00161 switch (Mt->give_dimension(eid))
00162 {
00163 case 1:
00164 l = Mt->give_length(eid);
00165 break;
00166 case 2:
00167 l = sqrt(Mt->give_area(eid));
00168 break;
00169 case 3:
00170 l = pow(Mt->give_volume(eid), 1.0/3.0);
00171 break;
00172 default:
00173 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00174 }
00175 tensor_vector_full(eps, epst, strain);
00176 mxv(dl,eps,sig);
00177 vector_tensor(sig,sigt,Mm->ip[ipp].ssst,stress);
00178
00179 for (i=0; i<3; i++)
00180 {
00181 xeq(i) = sigt(i,i)/dl(i,i) * l;
00182 if (xeq(i) < 0.0)
00183 xeq(i) = 0.0;
00184 }
00185
00186
00187
00188
00189
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 void fixortodam::compute_dam (vector &xeq, vector &omegao, vector &omega)
00204 {
00205 long i;
00206
00207 for (i=0; i<3; i++)
00208 {
00209 if (xeq(i) >= xeqtf(i))
00210 {
00211 omega(i) = 1.0;
00212 continue;
00213 }
00214
00215 if (xeq(i) <= xeqt0(i))
00216 {
00217 omega(i) = 0.0;
00218 continue;
00219 }
00220 omega(i)=xeqtf(i)*(xeq(i)-xeqt0(i))/(xeq(i)*(xeqtf(i)-xeqt0(i)));
00221 if (omega(i)<omegao(i))
00222 omega(i)=omegao(i);
00223 if (omega(i)<0)
00224 omega(i)=0;
00225 }
00226 }
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 void fixortodam::secstiffmat(matrix &d, long ipp, vector &omega)
00241 {
00242
00243
00244 Mm->loc_elmatstiff (d,ipp);
00245
00246 d(0,0)=pow(1-omega(0),2)*d(0,0);
00247 d(0,1)=(1-omega(0))*(1-omega(1))*d(0,1);
00248 d(0,2)=(1-omega(0))*(1-omega(2))*d(0,2);
00249 d(1,0)=d(0,1);
00250 d(1,1)=pow(1-omega(1),2)*d(1,1);
00251 d(1,2)=(1-omega(1))*(1-omega(2))*d(1,2);
00252 d(2,0)=d(0,2);
00253 d(2,1)=d(1,2);
00254 d(2,2)=pow(1-omega(2),2)*d(2,2);
00255 if (omega(0) == 1.0 && omega(1) == 1.0)
00256 d(3, 3) = 0;
00257 else
00258 d(3, 3) = pow((2 * (1 - omega(0))*(1 - omega(1))) / (2 - omega(0) - omega(1)), 2)*d(3,3);
00259 if (omega(0) == 1.0 && omega(2) == 1.0)
00260 d(4, 4) = 0;
00261 else
00262 d(4, 4) = pow((2 * (1 - omega(0))*(1 - omega(2))) / (2 - omega(0) - omega(2)), 2)*d(4,4);
00263 if (omega(1) == 1.0 && omega(2) == 1.0)
00264 d(5, 5) = 0;
00265 else
00266 d(5, 5) = pow((2 * (1 - omega(1))*(1 - omega(2))) / (2 - omega(1) - omega(2)), 2)*d(5,5);
00267 }
00268
00269
00270
00271 void fixortodam::nlstresses (long ipp, long im, long ido)
00272 {
00273 long i;
00274 long ncompstr=Mm->ip[ipp].ncompstr;
00275
00276 vector omega(3), omegao(3), xeq(3);
00277 vector epsg(ncompstr), eps(ncompstr), sig(ncompstr), sig_g(ncompstr);
00278 matrix epst(3,3), sigt(3,3), d(6,6);
00279 matrix tmat(3,3);
00280
00281
00282 for (i=0; i<ncompstr; i++)
00283 epsg(i) = Mm->ip[ipp].strain[i];
00284
00285
00286 Mm->loc_transf_mat(tmat, ipp);
00287
00288
00289 vector_tensor (epsg,epst,Mm->ip[ipp].ssst,strain);
00290 glmatrixtransf(epst, tmat);
00291 tensor_vector(eps,epst,Mm->ip[ipp].ssst,strain);
00292
00293
00294
00295
00296
00297 for (i=0; i<3; i++)
00298 omegao(i) = Mm->ip[ipp].eqother[ido+i];
00299
00300
00301 compute_eqdispl(ipp, epst, xeq);
00302
00303
00304 compute_dam(xeq, omegao, omega);
00305
00306
00307 secstiffmat(d, ipp, omega);
00308
00309
00310
00311 mxv(d, eps, sig);
00312
00313
00314 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00315 lgmatrixtransf(sigt, tmat);
00316 tensor_vector(sig_g,sigt,Mm->ip[ipp].ssst,stress);
00317
00318
00319 for (i=0; i<ncompstr; i++)
00320 Mm->ip[ipp].stress[i] = sig_g[i];
00321
00322
00323
00324 for (i=0; i<3; i++)
00325 Mm->ip[ipp].other[ido+i] = omega(i);
00326 for (i=0; i<3; i++)
00327 Mm->ip[ipp].other[ido+3+i] = xeq(i);
00328 }
00329
00330
00331
00332 void fixortodam::updateval (long ipp,long im,long ido)
00333 {
00334 long i;
00335
00336 for (i=0; i<6; i++)
00337 Mm->ip[ipp].eqother[ido+i] = Mm->ip[ipp].other[ido+i];
00338 }