00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include "layplate.h"
00005 #include "global.h"
00006 #include "vecttens.h"
00007 #include "intpoints.h"
00008 #include "matrix.h"
00009 #include "vector.h"
00010 #include "tensor.h"
00011 #include "elastisomat.h"
00012 #include "strretalg.h"
00013 #include "alias.h"
00014 #include "element.h"
00015
00016
00017
00018
00019 layplate::layplate (void)
00020 {
00021
00022 }
00023
00024
00025
00026
00027
00028 layplate::~layplate (void)
00029 {
00030
00031 }
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 void layplate::read (XFILE *in)
00043 {
00044 long i,j;
00045
00046
00047 xfscanf (in,"%k%ld", "num_lay", &nl);
00048
00049
00050
00051
00052 tm = new mattype*[nl];
00053 idm = new long*[nl];
00054 nm = new long[nl];
00055
00056 for (i=0;i<nl;i++){
00057 xfscanf (in,"%ld", &nm[i]);
00058 tm[i] = new mattype[nm[i]];
00059 idm[i] = new long[nm[i]];
00060
00061 for (j=0;j<nm[i];j++){
00062 xfscanf(in, "%k%m%k%ld", "mattype", &mattype_kwdset, (int*)(tm[i]+j), "num_inst", idm[i]+j);
00063 idm[i][j]--;
00064 }
00065 }
00066 xfscanf (in,"%k%ld", "num_iter", &nli);
00067 xfscanf (in,"%k%lf", "error", &err);
00068
00069 allocm(3,3,c);
00070 allocm(3,3,cn);
00071 }
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 void layplate::matstiff (matrix &d, long ipp,long ido)
00085 {
00086 long i,j,jj;
00087 double *th;
00088 double *zet;
00089 matrix e(3,3);
00090 matrix dn(3,3);
00091 matrix dnm(3,3);
00092 matrix dm(3,3);
00093 matrix pomat(3,3);
00094 matrix pommat(3,3);
00095 long eid = Mm->elip[ipp];
00096 double *k;
00097
00098 th = Mc->give_layer_thicke(eid);
00099 zet = Mc->give_layer_zcoord(eid);
00100
00101 fillm(0,e);
00102 fillm(0,d);
00103 backup(ipp,k,0);
00104 for (i=0;i<nl;i++){
00105 Mm->ip[ipp].tm = tm [i];
00106 Mm->ip[ipp].idm = idm[i];
00107 Mm->ip[ipp].nm = nm [i];
00108
00109 Mm->elmatstiff (e,ipp);
00110
00111 for (j=0;j<3;j++){
00112 for (jj=0;jj<3;jj++){
00113 dn[j][jj] += th[i]*e[j][jj];
00114 }
00115 }
00116
00117 for (j=0;j<3;j++){
00118 for (jj=0;jj<3;jj++){
00119 dnm[j][jj] += zet[i]*th[i]*e[j][jj];
00120 }
00121 }
00122
00123 for (j=0;j<3;j++){
00124 for (jj=0;jj<3;jj++){
00125 dm[j][jj] += th[i]*zet[i]*zet[i]*e[j][jj];
00126 }
00127 }
00128 }
00129 restore_values(ipp,k,0);
00130
00131 invm(dn,c,1.0e-20);
00132 cmulm(-1.0,dnm,pomat);
00133 mxm(c,pomat,cn);
00134 mxm(pomat,c,pommat);
00135 mxm(pommat,dnm,pomat);
00136 addm(pomat,dm,d);
00137
00138 }
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 void layplate::nlstresses (long ipp, long im, long ido)
00155 {
00156 long i,j,nc = Mm->ip[ipp].ncompstr;
00157 long eid = Mm->elip[ipp];
00158 double *k;
00159 double *layth;
00160 double *layz;
00161 double th;
00162 vector intfor(nc);
00163 vector df(3);
00164 vector eps(3);
00165 vector dpeps(3);
00166 vector deps(3);
00167 double *asd;
00168 double normf;
00169 FILE *out;
00170 char fname[14] = "kontrola.dat";
00171
00172 if (Mp->istep==97 && Mp->jstep==2 && Mm->elip[ipp]==3 && ipp==10)
00173 out = fopen(fname, "wt");
00174
00175 layth = Mc->give_layer_thicke(eid);
00176 layz = Mc->give_layer_zcoord(eid);
00177 for (i=0;i<3;i++)
00178 dpeps[i] = Mm->ip[ipp].eqother[i];
00179
00180 backup (ipp,k,1);
00181
00182 for (i=0;i<3;i++)
00183 for (j=0;j<3;j++)
00184 eps[i] += cn[i][j]*k[j];
00185
00186 for (i=0;i<3;i++)
00187 eps[i] += dpeps[i];
00188
00189 for (i=0;i<nli;i++)
00190 {
00191 stress_calc(df,eps,k,intfor,layz,layth,ipp,ido,out);
00192 if (Mp->istep==97 && Mp->jstep==2 && Mm->elip[ipp]==3 && ipp==10)
00193 fprintf(out," df= %lf\n",df[0]);
00194 normf = normv(df);
00195 if (normf<err)
00196 break;
00197 mxv(c,df,deps);
00198 for (j=0;j<3;j++)
00199 eps[j] -= deps[j];
00200 for (j=0;j<3;j++)
00201 Mm->ip[ipp].other[j] -= deps[j];
00202 if (i==nli-1 && normf>err)
00203 {
00204 print_err("cross-section equilibrium has not been reached (istep=%ld,jstep=%ld,nordf=%le, el=%ld, ipp=%ld)",
00205 __FILE__, __LINE__, __func__, Mp->istep, Mp->jstep, normf, Mm->elip[ipp], ipp);
00206 fprintf (Out,"\n cross-section equilibrium has not been reached");
00207 }
00208 }
00209 asd = Mm->ip[ipp].other;
00210 if (Mp->istep==97 && Mp->jstep==2 && Mm->elip[ipp]==3 && ipp==10)
00211 fclose(out);
00212
00213 Mm->storestress(0, ipp, intfor);
00214 restore_values (ipp,k,1);
00215 }
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 long layplate::compeqother(long ipp)
00229 {
00230 long sizeeqother=0;
00231 long i;
00232 mattype *btm = Mm->ip[ipp].tm;
00233 Mm->ip[ipp].ncompstr = 4;
00234
00235
00236 for (i=0;i<nl;i++){
00237 Mm->ip[ipp].tm = tm [i];
00238 sizeeqother += Mm->givencompeqother(ipp,0)+3;
00239 }
00240
00241 Mm->ip[ipp].tm = btm;
00242 Mm->ip[ipp].ncompstr = 3;
00243
00244 return sizeeqother+3;
00245 }
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 long layplate::compother(long ipp)
00259 {
00260 long sizeother=0;
00261 long i;
00262 mattype *btm = Mm->ip[ipp].tm;
00263 Mm->ip[ipp].ncompstr = 4;
00264
00265
00266 for (i=0;i<nl;i++){
00267 Mm->ip[ipp].tm = tm[i];
00268 sizeother += Mm->givencompother(ipp,0)+3;
00269 }
00270
00271 Mm->ip[ipp].tm = btm;
00272 Mm->ip[ipp].ncompstr = 3;
00273
00274 return sizeother+3;
00275 }
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 void layplate::backup(long ipp, double *&k, long j)
00286 {
00287 bcup.tm = Mm->ip[ipp].tm;
00288 bcup.idm = Mm->ip[ipp].idm;
00289 bcup.nm = Mm->ip[ipp].nm;
00290 bcup.ssst = Mm->ip[ipp].ssst;
00291 bcup.ncompstr = Mm->ip[ipp].ncompstr;
00292 Mm->ip[ipp].ssst = planestress;
00293 Mm->ip[ipp].ncompstr = 4;
00294
00295 if (j==1)
00296 {
00297 bcup.stress = Mm->ip[ipp].stress;
00298 k = Mm->ip[ipp].strain;
00299 Mm->ip[ipp].alloc_strains(1);
00300 Mm->ip[ipp].alloc_stresses(1);
00301 Mm->ip[ipp].clean_strains(1);
00302 Mm->ip[ipp].clean_stresses(1);
00303 }
00304 }
00305
00306
00307
00308
00309
00310
00311
00312
00313 void layplate::restore_values(long ipp, double *k, long j)
00314 {
00315 long i;
00316
00317 Mm->ip[ipp].tm = bcup.tm;
00318 Mm->ip[ipp].idm = bcup.idm;
00319 Mm->ip[ipp].nm = bcup.nm;
00320 Mm->ip[ipp].ssst = bcup.ssst;
00321 Mm->ip[ipp].ncompstr = bcup.ncompstr;
00322
00323 if (j==1)
00324 {
00325 delete [] Mm->ip[ipp].strain;
00326 Mm->ip[ipp].strain = k;
00327 for (i=0; i<3; i++)
00328 bcup.stress[i] = Mm->ip[ipp].stress[i];
00329 delete [] Mm->ip[ipp].stress;
00330 Mm->ip[ipp].stress = bcup.stress;
00331 }
00332 }
00333
00334
00335
00336
00337 void layplate::stress_calc(vector &df, vector &eps, double *k, vector &intfor, double *layz, double *layth, long ipp, long ido, FILE *&out)
00338 {
00339 long i,j,nc = Mm->ip[ipp].ncompstr;
00340 long neq;
00341
00342 for (i=0;i<3;i++){
00343 df[i] = 0.0;
00344 }
00345
00346 ido = 3;
00347
00348 for (i=0;i<nl;i++)
00349 {
00350 Mm->ip[ipp].tm = tm [i];
00351 Mm->ip[ipp].idm = idm[i];
00352 Mm->ip[ipp].nm = nm [i];
00353
00354 for (j=0;j<3;j++){
00355 Mm->ip[ipp].strain[j] = k[j]*layz[i]+eps[j];
00356 }
00357
00358 Mm->computenlstresses(ipp,0,ido);
00359
00360 for (j=0;j<3;j++){
00361 intfor[j] += Mm->ip[ipp].stress[j]*layth[i]*layz[i];
00362 }
00363
00364 neq = Mm->givencompeqother(ipp,0);
00365 for (j=0;j<3;j++){
00366 Mm->ip[ipp].other[ido+neq+j] = Mm->ip[ipp].stress[j];
00367 }
00368
00369 if (Mp->istep==97 && Mp->jstep==2 && Mm->elip[ipp]==3 && ipp==10)
00370 fprintf(out," %lf %lf",Mm->ip[ipp].stress[0],Mm->ip[ipp].strain[0]);
00371
00372 for (j=0;j<3;j++)
00373 df[j] += Mm->ip[ipp].stress[j]*layth[i];
00374
00375 ido += neq+3;
00376 }
00377 }
00378
00379
00380
00381
00382
00383
00384 void layplate::updateval (long ipp, long im, long ido)
00385 {
00386 long i,n = Mm->givencompeqother(ipp,im);
00387
00388 for (i=0;i<n;i++)
00389 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00390 }
00391