00001 #include "homog.h"
00002 #include "global.h"
00003 #include "element.h"
00004 #include "intpoints.h"
00005 #include "globmat.h"
00006 #include "gtopology.h"
00007 #include "math.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 int C012(int x){
00019 if (x<0){
00020 x+=3;
00021 }
00022 else if(x>2){
00023 x-=3;
00024 }
00025 return x;
00026 }
00027
00028 void homogenization (FILE *out, long lcid)
00029 {
00030
00031 if (Mp->homog!=3){
00032
00033
00034
00035
00036
00037
00038
00039 long i, j, k, ncomp, id, ipp, tnipe, ndofn, counter;
00040 double *stress_av, *strain_av;
00041 double *stress_max, *strain_max;
00042 double *r;
00043 double x_sum, y_sum, z_max, z_sum;
00044 double stress_eq_max=0;
00045 double mu[3], E[9], nu[9];
00046 double G_en, nu_en, E_en, W_en, help_a, help_b;
00047 FILE *propfile;
00048 char filename[]="Elas_1.out";
00049
00050 ncomp=6;
00051
00052 stress_av = new double [ncomp];
00053 strain_av = new double [ncomp];
00054
00055 stress_max = new double [ncomp];
00056 strain_max = new double [ncomp];
00057
00058 if((propfile=fopen(filename, "a"))==NULL){
00059 printf("error in opening %s file for writing\n", filename);
00060 perror("");
00061 }
00062
00063 for(i=0;i<ncomp;i++){
00064 stress_av[i]=0;
00065 strain_av[i]=0;
00066 stress_max[i]=0;
00067 strain_max[i]=0;
00068 }
00069
00070
00071
00072 for (i=0; i<Mt->ne; i++){
00073 ipp = Mt->elements[i].ipp[0][0];
00074 tnipe = Mt->give_totnip(i);
00075 ncomp = Mm->ip[ipp].ncompstr;
00076 id = lcid*ncomp;
00077 for (j=0; j<tnipe; j++){
00078 for (k=0; k<ncomp; k++){
00079
00080 help_a=Mm->ip[ipp+j].stress[id+k];
00081 help_b=Mm->ip[ipp+j].strain[id+k]-Mm->eigstrains[j][k];
00082
00083 stress_av[k]+=help_a;
00084 strain_av[k]+=help_b;
00085
00086
00087
00088
00089 if(help_a<stress_max[k]){
00090 stress_max[k]=help_a;
00091 }
00092 if(help_b<strain_max[k]){
00093 strain_max[k]=help_b;
00094 }
00095 }
00096
00097 help_a=sqrt(2.)/2.*
00098 sqrt((Mm->ip[ipp+j].stress[id+0]-Mm->ip[ipp+j].stress[id+1])*(Mm->ip[ipp+j].stress[id+0]-Mm->ip[ipp+j].stress[id+1])+
00099 (Mm->ip[ipp+j].stress[id+1]-Mm->ip[ipp+j].stress[id+2])*(Mm->ip[ipp+j].stress[id+1]-Mm->ip[ipp+j].stress[id+2])+
00100 (Mm->ip[ipp+j].stress[id+2]-Mm->ip[ipp+j].stress[id+0])*(Mm->ip[ipp+j].stress[id+2]-Mm->ip[ipp+j].stress[id+0])+
00101 6*Mm->ip[ipp+j].stress[id+3]*Mm->ip[ipp+j].stress[id+3]+
00102 6*Mm->ip[ipp+j].stress[id+4]*Mm->ip[ipp+j].stress[id+4]+
00103 6*Mm->ip[ipp+j].stress[id+5]*Mm->ip[ipp+j].stress[id+5]);
00104
00105 if(help_a>stress_eq_max){
00106 stress_eq_max=help_a;
00107 }
00108
00109 }
00110 }
00111
00112
00113 for(i=0;i<ncomp;i++){
00114
00115 strain_av[i]/=(-Mt->ne*tnipe);
00116 stress_av[i]/=(-Mt->ne*tnipe);
00117 }
00118
00119
00120 if (Mp->homog==1){
00121
00122
00123 G_en=(strain_av[3]*stress_av[3]+strain_av[4]*stress_av[4]+
00124 strain_av[5]*stress_av[5])/
00125 (strain_av[3]*strain_av[3]+strain_av[4]*strain_av[4]+
00126 strain_av[5]*strain_av[5]);
00127
00128 help_a=(strain_av[0]*strain_av[0]+strain_av[1]*strain_av[1]+
00129 strain_av[2]*strain_av[2]);
00130
00131 help_b=2*(strain_av[0]*strain_av[1]+strain_av[0]*strain_av[2]+
00132 strain_av[1]*strain_av[2]);
00133
00134
00135 W_en=0.5*(strain_av[0]*stress_av[0]+strain_av[1]*stress_av[1]+
00136 strain_av[2]*stress_av[2]);
00137
00138 nu_en=(W_en-G_en*help_a)/(G_en*(help_b-help_a)+2*W_en);
00139 E_en=2*G_en*(1+nu_en);
00140
00141
00142 fprintf(propfile, "ENER. %f %f %f %f %f ", E_en, nu_en, G_en, E_en/(3.*(1.-2*nu_en)), W_en);
00143
00144 fprintf(propfile, "STRE ");
00145 for(i=0;i<ncomp;i++){
00146 fprintf(propfile,"%f ", stress_av[i]);
00147 }
00148
00149 fprintf(propfile,"STRN ");
00150 for(i=0;i<ncomp;i++){
00151 fprintf(propfile,"%f ", strain_av[i]);
00152 }
00153
00154
00155 fprintf(propfile,"MHH ");
00156 fprintf(propfile,"%f ", stress_eq_max);
00157
00158
00159 fprintf(propfile,"MAX_STRE ");
00160 for(i=0;i<ncomp;i++){
00161 fprintf(propfile,"%f ", stress_max[i]);
00162 }
00163
00164 fprintf(propfile,"MAX_STRN ");
00165 for(i=0;i<ncomp;i++){
00166 fprintf(propfile,"%f ", strain_max[i]);
00167 }
00168
00169
00170 mu[0]=stress_av[3]/strain_av[3];
00171 mu[1]=stress_av[4]/strain_av[4];
00172 mu[2]=stress_av[5]/strain_av[5];
00173
00174
00175 fprintf(propfile, "E ");
00176 for (i=0; i<=2; i++){
00177 for (j=0; j<=2; j++){
00178 E[3*i+j]=(stress_av[C012(i+0)]+stress_av[C012(i+1)]+stress_av[C012(i+2)])/
00179 (strain_av[C012(i+0)]+(stress_av[C012(i+1)]+stress_av[C012(i+2)])/(2*mu[j]));
00180 fprintf(propfile, "%f ", E[3*i+j]);
00181 }
00182 }
00183
00184
00185 fprintf(propfile, "Nu ");
00186 for (i=0; i<=2; i++){
00187 for (j=0; j<=2; j++){
00188 nu[3*i+j]=(E[3*i+j]-2*mu[i])/(2*mu[i]);
00189 fprintf(propfile, "%f ", nu[3*i+j]);
00190 }
00191 }
00192
00193 fprintf(propfile, "\n");
00194
00195 }
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 else if(Mp->homog==2){
00248
00249
00250 z_max=0.;
00251 x_sum=0.;
00252 y_sum=0.;
00253 z_sum=0.;
00254 for (i=0; i<Mt->nn; i++){
00255 if(Gtm->gnodes[i].z>z_max){
00256 z_max=Gtm->gnodes[i].z;
00257 }
00258 }
00259
00260 printf("\n ");
00261
00262 for(i=0; i<Mt->nn; i++){
00263 if(z_max==Gtm->gnodes[i].z){
00264 ndofn = Mt->give_ndofn(i);
00265 r = new double[ndofn];
00266 noddispl (lcid, r, i);
00267 x_sum+=r[0];
00268 y_sum+=r[1];
00269 z_sum+=r[2];
00270
00271 delete [] r;
00272 }
00273 }
00274
00275 x_sum/=(z_max*z_max);
00276 y_sum/=(z_max*z_max);
00277 z_sum/=(z_max*z_max);
00278
00279 fprintf(propfile, "z_strn %f x_strn %f y_strn %f av_displ_z %f ", z_sum/(z_max-1), x_sum/(z_max-1), y_sum/(z_max-1), z_sum);
00280
00281 for(i=0;i<ncomp;i++){
00282 stress_av[i]=0;
00283 strain_av[i]=0;
00284 stress_max[i]=0;
00285 strain_max[i]=0;
00286 }
00287
00288
00289
00290
00291 for (k=0; k<ncomp; k++){
00292 stress_av[k]=0;
00293 strain_av[k]=0;
00294 }
00295 counter=0;
00296
00297 for (i=0; i<Mt->ne; i++){
00298
00299 if (Mt->elements[i].idm[0]==(28-1)){
00300 ipp = Mt->elements[i].ipp[0][0];
00301 tnipe = Mt->give_totnip(i);
00302 ncomp = Mm->ip[ipp].ncompstr;
00303 id = lcid*ncomp;
00304 for (j=0; j<tnipe; j++){
00305 counter++;
00306 for (k=0; k<ncomp; k++){
00307
00308 help_a=Mm->ip[ipp+j].stress[id+k];
00309 help_b=Mm->ip[ipp+j].strain[id+k];
00310
00311 stress_av[k]+=help_a;
00312 strain_av[k]+=help_b;
00313 }
00314 }
00315 }
00316 }
00317
00318
00319 for(i=0;i<ncomp;i++){
00320
00321 strain_av[i]/=(counter);
00322 stress_av[i]/=(counter);
00323 }
00324
00325 fprintf(propfile,"SIG_POR ");
00326
00327
00328 for(i=0;i<ncomp;i++){
00329 fprintf(propfile,"%f ", stress_av[i]);
00330 }
00331
00332 fprintf(propfile,"EPS_POR ");
00333
00334 for(i=0;i<ncomp;i++){
00335 fprintf(propfile,"%f ", strain_av[i]);
00336 }
00337
00338
00339 for (k=0; k<ncomp; k++){
00340 stress_av[k]=0;
00341 strain_av[k]=0;
00342 }
00343 counter=0;
00344
00345 for (i=0; i<Mt->ne; i++){
00346
00347 if(Mt->elements[i].idm[0]!=(28-1) && Mt->elements[i].idm[0]!=(26-1) && Mt->elements[i].idm[0]!=(29-1)){
00348 ipp = Mt->elements[i].ipp[0][0];
00349 tnipe = Mt->give_totnip(i);
00350 ncomp = Mm->ip[ipp].ncompstr;
00351 id = lcid*ncomp;
00352 for (j=0; j<tnipe; j++){
00353 counter++;
00354 for (k=0; k<ncomp; k++){
00355
00356 help_a=Mm->ip[ipp+j].stress[id+k];
00357 help_b=Mm->ip[ipp+j].strain[id+k];
00358 stress_av[k]+=help_a;
00359 strain_av[k]+=help_b;
00360 }
00361 }
00362 }
00363 }
00364
00365
00366 fprintf(propfile,"SIG_SOL_RVE ");
00367
00368 for(i=0;i<ncomp;i++){
00369 fprintf(propfile,"%f ", stress_av[i]/(z_max-1.)/(z_max-1.)/(z_max-1.)/8.);
00370 }
00371
00372 fprintf(propfile,"EPS_SOL_RVE ");
00373
00374 for(i=0;i<ncomp;i++){
00375 fprintf(propfile,"%f ", strain_av[i]/(z_max-1.)/(z_max-1.)/(z_max-1.)/8.);
00376 }
00377
00378 fprintf(propfile,"SIG_SOL ");
00379
00380 for(i=0;i<ncomp;i++){
00381 fprintf(propfile,"%f ", stress_av[i]/counter);
00382 }
00383
00384 fprintf(propfile,"EPS_SOL ");
00385
00386 for(i=0;i<ncomp;i++){
00387 fprintf(propfile,"%f ", strain_av[i]/counter);
00388 }
00389
00390
00391 fprintf(propfile,"\n");
00392
00393 }
00394
00395 fclose(propfile);
00396
00397 delete [] stress_av;
00398 delete [] strain_av;
00399 delete [] stress_max;
00400 delete [] strain_max;
00401
00402 }
00403 }