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 vector intfor(nc);
00162 vector df(3);
00163 vector eps(3);
00164 vector dpeps(3);
00165 vector deps(3);
00166 double normf;
00167
00168
00169
00170 Mm->ip[ipp].other[3]=1.0;
00171 layth = Mc->give_layer_thicke(eid);
00172 layz = Mc->give_layer_zcoord(eid);
00173
00174
00175
00176 backup (ipp,k,1);
00177
00178 for (i=0;i<3;i++)
00179 for (j=0;j<3;j++)
00180 eps[i] += cn[i][j]*k[j];
00181
00182
00183
00184
00185 for (i=0;i<nli;i++)
00186 {
00187 stress_calc(df,eps,k,intfor,layz,layth,ipp,ido);
00188 normf = normv(df);
00189 if (normf<err)
00190 break;
00191 mxv(c,df,deps);
00192 for (j=0;j<3;j++)
00193 eps[j] -= deps[j];
00194
00195
00196 if (i==nli-1 && normf>err)
00197 {
00198 Mm->ip[ipp].other[3]=0.5;
00199 }
00200
00201
00202 }
00203
00204 Mm->storestress(0, ipp, intfor);
00205 restore_values (ipp,k,1);
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219 long layplate::compeqother(long ipp)
00220 {
00221 long sizeeqother=0;
00222 long i;
00223 mattype *btm = Mm->ip[ipp].tm;
00224 Mm->ip[ipp].ncompstr = 4;
00225
00226
00227 for (i=0;i<nl;i++){
00228 Mm->ip[ipp].tm = tm [i];
00229 sizeeqother += Mm->givencompeqother(ipp,0)+3;
00230 }
00231
00232 Mm->ip[ipp].tm = btm;
00233 Mm->ip[ipp].ncompstr = 3;
00234
00235 return sizeeqother+5;
00236 }
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 long layplate::compother(long ipp)
00250 {
00251 long sizeother=0;
00252 long i;
00253 mattype *btm = Mm->ip[ipp].tm;
00254 Mm->ip[ipp].ncompstr = 4;
00255
00256
00257 for (i=0;i<nl;i++){
00258 Mm->ip[ipp].tm = tm[i];
00259 sizeother += Mm->givencompother(ipp,0)+3;
00260 }
00261
00262 Mm->ip[ipp].tm = btm;
00263 Mm->ip[ipp].ncompstr = 3;
00264
00265 return sizeother+5;
00266 }
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276 void layplate::backup(long ipp, double *&k, long j)
00277 {
00278 bcup.tm = Mm->ip[ipp].tm;
00279 bcup.idm = Mm->ip[ipp].idm;
00280 bcup.nm = Mm->ip[ipp].nm;
00281 bcup.ssst = Mm->ip[ipp].ssst;
00282 bcup.ncompstr = Mm->ip[ipp].ncompstr;
00283 Mm->ip[ipp].ssst = planestress;
00284 Mm->ip[ipp].ncompstr = 4;
00285
00286 if (j==1)
00287 {
00288 bcup.stress = Mm->ip[ipp].stress;
00289 k = Mm->ip[ipp].strain;
00290 Mm->ip[ipp].alloc_strains(1);
00291 Mm->ip[ipp].alloc_stresses(1);
00292 Mm->ip[ipp].clean_strains(1);
00293 Mm->ip[ipp].clean_stresses(1);
00294 }
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304 void layplate::restore_values(long ipp, double *k, long j)
00305 {
00306 long i;
00307
00308 Mm->ip[ipp].tm = bcup.tm;
00309 Mm->ip[ipp].idm = bcup.idm;
00310 Mm->ip[ipp].nm = bcup.nm;
00311 Mm->ip[ipp].ssst = bcup.ssst;
00312 Mm->ip[ipp].ncompstr = bcup.ncompstr;
00313 bcup.tm = NULL;
00314 bcup.idm = NULL;
00315
00316 if (j==1)
00317 {
00318 delete [] Mm->ip[ipp].strain;
00319 Mm->ip[ipp].strain = k;
00320 for (i=0; i<3; i++)
00321 bcup.stress[i] = Mm->ip[ipp].stress[i];
00322 delete [] Mm->ip[ipp].stress;
00323 Mm->ip[ipp].stress = bcup.stress;
00324 bcup.stress = NULL;
00325 }
00326 }
00327
00328
00329
00330
00331 void layplate::stress_calc(vector &df, vector &eps, double *k, vector &intfor, double *layz, double *layth, long ipp, long ido)
00332 {
00333 long i,j;
00334 long neq;
00335 double *stra;
00336
00337 for (i=0;i<3;i++){
00338 df[i] = 0.0;
00339 intfor[i] = 0.0;
00340 }
00341
00342 ido = 5;
00343 Mm->ip[ipp].other[4]=0.0;
00344
00345 for (i=0;i<nl;i++)
00346 {
00347 Mm->ip[ipp].tm = tm [i];
00348 Mm->ip[ipp].idm = idm[i];
00349 Mm->ip[ipp].nm = nm [i];
00350
00351 for (j=0;j<3;j++){
00352 Mm->ip[ipp].strain[j] = k[j]*layz[i]+eps[j];
00353 }
00354 stra = Mm->ip[ipp].strain;
00355
00356 Mm->computenlstresses(ipp,0,ido);
00357
00358 for (j=0;j<3;j++){
00359 intfor[j] += Mm->ip[ipp].stress[j]*layth[i]*layz[i];
00360 }
00361
00362 neq = Mm->givencompeqother(ipp,0);
00363 for (j=0;j<3;j++){
00364 Mm->ip[ipp].other[ido+neq+j] = Mm->ip[ipp].stress[j];
00365 }
00366
00367 for (j=0;j<3;j++)
00368 df[j] += Mm->ip[ipp].stress[j]*layth[i];
00369
00370 for (j=0;j<3;j++)
00371 Mm->ip[ipp].other[4] += Mm->ip[ipp].other[ido+4];
00372
00373 ido += neq+3;
00374 }
00375 }
00376
00377
00378
00379
00380
00381
00382 void layplate::updateval (long ipp, long im, long ido)
00383 {
00384 long i,n = Mm->givencompeqother(ipp,im);
00385
00386 for (i=0;i<n;i++)
00387 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00388 }
00389
00390 double layplate::dstep_red(long ipp, long im, long ido)
00391 {
00392 return (Mm->ip[ipp].other[3]);
00393 }
00394