00001 #include "soilbeam.h"
00002 #include "global.h"
00003 #include "globmat.h"
00004 #include "genfile.h"
00005 #include "node.h"
00006 #include "element.h"
00007 #include "beamel3d.h"
00008 #include "intpoints.h"
00009 #include "elastisomat.h"
00010 #include <stdlib.h>
00011 #include <math.h>
00012
00013 soilbeam::soilbeam (void)
00014 {
00015 long i;
00016
00017
00018 nne=2;
00019
00020 ndofe=12;
00021
00022 tncomp=6;
00023
00024 napfun=3;
00025
00026 intordmm=4;
00027
00028 intordism=2;
00029
00030 ssst=spacebeam;
00031
00032
00033 nb=1;
00034
00035
00036 ncomp = new long [nb];
00037 ncomp[0]=1;
00038
00039
00040 cncomp = new long [nb];
00041 cncomp[0]=0;
00042
00043
00044
00045 nip = new long* [nb];
00046 intordsm = new long* [nb];
00047 for (i=0;i<nb;i++){
00048 nip[i] = new long [nb];
00049 intordsm[i] = new long [nb];
00050 }
00051 nip[0][0]=2;
00052 intordsm[0][0]=2;
00053
00054 c1= new double [3];
00055 c2= new double [3];
00056
00057
00058 bPod=1.;
00059 tnip=2;
00060
00061 }
00062
00063 soilbeam::~soilbeam (void)
00064 {
00065 long i;
00066
00067 for (i=0;i<nb;i++){
00068 delete [] nip[i];
00069 }
00070 delete [] nip;
00071 }
00072
00073
00074 void soilbeam::transf_matrix (ivector &nodes,matrix &tmat)
00075 {
00076 long i,i6,n,m;
00077
00078 fillm (0.0,tmat);
00079
00080 n=nodes.n;
00081 m=tmat.m;
00082 for (i=0;i<m;i++){
00083 tmat[i][i]=1.0;
00084 }
00085
00086 for (i=0;i<n;i++){
00087 if (Mt->nodes[nodes[i]].transf>0){
00088 i6=i*6;
00089 tmat[i6][i6] = Mt->nodes[nodes[i]].e1[0]; tmat[i6][i6+1] = Mt->nodes[nodes[i]].e2[0]; tmat[i6][i6+2] = Mt->nodes[nodes[i]].e3[0];
00090 tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1]; tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1]; tmat[i6+1][i6+2] = Mt->nodes[nodes[i]].e3[1];
00091 tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2]; tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2]; tmat[i6+2][i6+2] = Mt->nodes[nodes[i]].e3[2];
00092 i6=i*6+3;
00093 tmat[i6][i6] = Mt->nodes[nodes[i]].e1[0]; tmat[i6][i6+1] = Mt->nodes[nodes[i]].e2[0]; tmat[i6][i6+2] = Mt->nodes[nodes[i]].e3[0];
00094 tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1]; tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1]; tmat[i6+1][i6+2] = Mt->nodes[nodes[i]].e3[1];
00095 tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2]; tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2]; tmat[i6+2][i6+2] = Mt->nodes[nodes[i]].e3[2];
00096
00097 }
00098 }
00099 }
00100 void soilbeam::beam_transf_matrix (matrix &tmat,double &dl,vector &vec,vector &x,vector &y,vector &z,long eid)
00101 {
00102 double c;
00103 fillm (0.0,tmat);
00104
00105
00106
00107
00108
00109
00110
00111
00112 tmat[0][0]=x[1]-x[0];
00113 tmat[1][0]=y[1]-y[0];
00114 tmat[2][0]=z[1]-z[0];
00115 dl=sqrt((tmat[0][0]*tmat[0][0]+tmat[1][0]*tmat[1][0]+tmat[2][0]*tmat[2][0]));
00116 if (dl<Mp->zero){
00117 fprintf (stderr,"\n\n zero length of the %ld beamel3d element",eid);
00118 fprintf (stderr,"\n in function beamel3d::beam_transf_matrix.\n");
00119 }
00120 tmat[0][0]=tmat[0][0]/dl;
00121 tmat[1][0]=tmat[1][0]/dl;
00122 tmat[2][0]=tmat[2][0]/dl;
00123
00124 tmat[0][1] =vec[1]*tmat[2][0]-vec[2]*tmat[1][0];
00125 tmat[1][1] =vec[2]*tmat[0][0]-vec[0]*tmat[2][0];
00126 tmat[2][1] =vec[0]*tmat[1][0]-vec[1]*tmat[0][0];
00127 c=sqrt((tmat[0][1]*tmat[0][1]+tmat[1][1]*tmat[1][1]+tmat[2][1]*tmat[2][1]));
00128 if (c<Mp->zero){
00129 fprintf (stderr,"\n\n zero vec of the %ld beamel3d element",eid);
00130 fprintf (stderr,"\n in function beamel3d::beam_transf_matrix.\n");
00131 }
00132 tmat[0][1]=tmat[0][1]/c;
00133 tmat[1][1]=tmat[1][1]/c;
00134 tmat[2][1]=tmat[2][1]/c;
00135
00136 tmat[0][2]=tmat[1][0]*tmat[2][1]-tmat[2][0]*tmat[1][1];
00137 tmat[1][2]=tmat[2][0]*tmat[0][1]-tmat[0][0]*tmat[2][1];
00138 tmat[2][2]=tmat[0][0]*tmat[1][1]-tmat[1][0]*tmat[0][1];
00139 c=sqrt((tmat[0][2]*tmat[0][2]+tmat[1][2]*tmat[1][2]+tmat[2][2]*tmat[2][2]));
00140 if (c<Mp->zero){
00141 fprintf (stderr,"\n\n zero vec of the %ld beamel3d element",eid);
00142 fprintf (stderr,"\n in function beamel2d::beam_transf_matrix.\n");
00143 }
00144 tmat[0][2]=tmat[0][2]/c;
00145 tmat[1][2]=tmat[1][2]/c;
00146 tmat[2][2]=tmat[2][2]/c;
00147 }
00148 void soilbeam::geom_matrix (matrix &n,double s,double dl,double gy,double gz)
00149 {
00150 double ll,aj1;
00151 ll=dl*dl;
00152
00153 n[0][0] = 1.-s;
00154 n[0][6] = s;
00155
00156 aj1=1./(1.+2.*gy);
00157 n[2][2] =aj1*(1.+2.*gy-2.*gy*s-3.*s*s+2*s*s*s);
00158 n[2][4] =aj1*(-(1.+gy)*s+(2.+gy)*s*s-s*s*s)*dl;
00159 n[2][8] =aj1*(2.*gy*s+3.*s*s-2.*s*s*s);
00160 n[2][10]=aj1*(gy*s+(1.-gy)*s*s-s*s*s)*dl;
00161
00162 n[4][2] = aj1*(6.*s-6.*s*s)/dl;
00163 n[4][4] = aj1*(1.+2.*gy-2.*(2.+gy)*s+3.*s*s);
00164 n[4][8] =-n[4][2];
00165 n[4][10]= aj1*(-2.*(1.-gy)*s+3.*s*s);
00166
00167
00168 aj1=1./(1.+2.*gz);
00169 n[1][1] = aj1*(1.+2.*gz-2.*gz*s-3.*s*s+2.*s*s*s);
00170 n[1][5] =-aj1*(-(1.+gz)*s+(2.+gz)*s*s-s*s*s)*dl;
00171 n[1][7] = aj1*(2.*gz*s+3.*s*s-2.*s*s*s);
00172 n[1][11]=-aj1*(gz*s+(1.-gz)*s*s-s*s*s)*dl;
00173
00174 n[5][1] =-aj1*(6.*s-6.*s*s)/dl;
00175 n[5][5] = aj1*(1.+2.*gz-2.*(2.+gz)*s+3.*s*s);
00176 n[5][7] =-n[5][1];
00177 n[5][11]= aj1*(-2.*(1.-gz)*s+3.*s*s);
00178
00179 n[3][3] = 1.-s;
00180 n[3][9] = s;
00181
00182 }
00183
00184
00185 void soilbeam::res_stiffness_matrix (long eid,matrix &sm)
00186 {
00187 stiffness_matrix (eid,0,0,sm);
00188 }
00189
00190 void soilbeam::stiffness_matrix (long eid,long ri,long ci,matrix &sm)
00191 {
00192 long i,ipp,i1,transf, imat;
00193 double a2,a4,e,g,a,ixyz[3],beta[2],dl,ll,gy,gz,g2;
00194 ivector nodes(nne);
00195 vector vec(3),x(nne),y(nne),z(nne);
00196 matrix c(tncomp,tncomp),tran(3,3);
00197 Mt->give_elemnodes (eid,nodes);
00198 Mt->give_node_coord3d (x,y,z,eid);
00199 Mc->give_vectorlcs (eid, vec);
00200 beam_transf_matrix (tran,dl, vec,x,y,z,eid);
00201
00202 fillm (0.0,sm);
00203 ipp=Mt->elements[eid].ipp[ri][ci];
00204 Mm->matstiff (c,ipp);
00205 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00206 Mc->give_area (eid,a);
00207 Mc->give_mominer (eid,ixyz);
00208 Mc->give_shearcoeff (eid,beta);
00209
00210 e=Mm->eliso[imat].e; g=e/(1.+2.*Mm->eliso[imat].nu);
00211
00212
00213 c1[0]=c[0][0];
00214 c1[1]=c[1][1];
00215 c1[2]=c[2][2];
00216 c2[0]=c[3][3];
00217 c2[1]=c[4][4];
00218 c2[2]=c[5][5];
00219 ll=dl*dl;
00220
00221 if(beta[0]<Mp->zero) gy=0.0;
00222 else gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
00223
00224 if(beta[1]<Mp->zero) gz=0.;
00225 else gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
00226
00227
00228 sm[0][0]= c1[0]*dl/3;
00229 sm[0][6]= c1[0]*dl/6;
00230 g2=(c2[1]+c2[2])/2.;
00231 a2=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt((c1[1]+c1[2])/2.*g2) );
00232 a4=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt(g2*g2*g2/((c1[1]+c1[2])/2.)) );
00233 sm[3][3]= a2/3.*dl+a4/dl;
00234 sm[3][9]= a2/6.*dl-a4/dl;
00235
00236
00237 g2=gy*gy;
00238 a2=bPod*(c1[2]+sqrt(c1[2]*c2[2])*2./bPod) *dl/(1.+2.*gy)/(1.+2.*gy);
00239 a4=bPod*(c2[2]+sqrt(c2[2]*c2[2]*c2[2]/c1[2])/bPod);
00240
00241 sm[2][2] = a2*(4.*g2/3.+1.4*gy+13./35.)+a4*4.*(g2+gy+.3);
00242 sm[2][4] =-(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl;
00243 sm[2][8] = a2*(2./3.*g2+0.6*gy+9./70.)-a4*4*(g2+gy+.3);
00244 sm[2][10]= (a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl;
00245 sm[4][4] = (a2*(g2/30.+gy/30.+1./105.)+a4*(g2+gy+0.4)/3.)*ll;
00246 sm[4][8] =-(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl;
00247 sm[4][10]=-(a2*(g2/30.+gy/30.+1./140.)+a4*(g2+gy+0.1)/3.)*ll;
00248
00249 g2=gz*gz;
00250 a2=bPod*(c1[1]+sqrt(c1[1]*c2[1])*2./bPod)*dl/(1.+2.*gz)/(1.+2.*gz);
00251 a4=bPod*(c2[1]+sqrt(c2[1]*c2[1]*c2[1]/c1[1])/bPod);
00252
00253 sm[1][1] = a2*(4.*g2/3.+1.4*gz+13./35.)+a4*4.*(g2+gz+.3);
00254 sm[1][5] = (a2*(g2/6.+11.*gz/60.+11./210.)+a4/10. )*dl;
00255 sm[1][7] = a2*(2./3.*g2+0.6*gz+9./70.)-a4*4.*(g2+gz+.3);
00256 sm[1][11]=-(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl;
00257 sm[5][5] = (a2*(g2/30.+gz/30.+1./105.)+a4*(g2+gz+0.4)/3.)*ll;
00258 sm[5][11]=-(a2*(g2/30.+gz/30+1./140.)+a4*(g2+gz+0.1)/3.)*ll;
00259 sm[5][7] = (a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl;
00260
00261 for (i=0;i<6;i++){
00262 i1=i+6;
00263 sm[i1][i1]=sm[i][i];
00264 sm[i1][i]=sm[i][i1];
00265 }
00266
00267 sm[4][2] = sm[2][4];
00268 sm[5][1] = sm[1][5];
00269 sm[7][5] = sm[5][7];
00270 sm[8][4] = sm[4][8];
00271 sm[11][1]= sm[1][11];
00272 sm[10][2]= sm[2][10];
00273 sm[7][11]=-sm[1][5];
00274 sm[11][7]=-sm[1][5];
00275 sm[8][10]=-sm[2][4];
00276 sm[10][8]=-sm[2][4];
00277
00278
00279 lgmatrixtransfblock (sm,tran);
00280
00281 transf = Mt->locsystems (nodes);
00282 if (transf>0){
00283 matrix tmat (ndofe,ndofe);
00284 transf_matrix (nodes,tmat);
00285 glmatrixtransf (sm,tmat);
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 }
00298
00299
00300
00301 void soilbeam::strains (long lcid,long eid,long ri,long ci)
00302 {
00303 long ipp,i,j,ij,transf;
00304 double dl;
00305 ivector nodes(nne);
00306 vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),eps(6);
00307 matrix tran(3,3);
00308
00309
00310
00311 ipp=Mt->elements[eid].ipp[ri][ci];
00312 Mt->give_elemnodes (eid,nodes);
00313 Mt->give_node_coord3d (x,y,z,eid);
00314 eldispl (lcid,eid,rl.a);
00315 Mc->give_vectorlcs (eid, vec);
00316 beam_transf_matrix (tran,dl, vec,x,y,z,eid);
00317
00318
00319 transf = Mt->locsystems (nodes);
00320 if (transf>0){
00321 matrix tmat (ndofe,ndofe);
00322 transf_matrix (nodes,tmat);
00323
00324 lgvectortransf (rg,rl,tmat);
00325 }
00326 else{
00327 copyv (rl,rg);
00328 }
00329
00330
00331
00332 ij=0;
00333 for (i=0;i<intordsm[0][0];i++){
00334 for (j=0;j<6;j++){
00335 eps[j]=rl[ij];
00336 ij++;
00337 }
00338
00339
00340 Mm->storestrain(lcid,ipp,eps);
00341 ipp++;
00342 }
00343 }
00344
00345 void soilbeam::res_internal_forces (long lcid,long eid,vector &ifor)
00346 {
00347 internal_forces (lcid, eid, 0, 0, ifor);
00348 }
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361 void soilbeam::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
00362 {
00363 long ipp,imat,i,j,k;
00364 double dl,ll,a2,a4,ixyz[3],beta[2],a,e,g,gy,gz,g2;
00365 ivector nodes(nne);
00366 vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),f(ndofe),eps(6);
00367 matrix d(tncomp,tncomp),c(tncomp,tncomp),tran(3,3),r(2,ndofe);
00368
00369 Mt->give_elemnodes (eid,nodes);
00370 Mt->give_node_coord3d (x,y,z,eid);
00371
00372 ipp=Mt->elements[eid].ipp[ri][ci];
00373 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00374 Mm->matstiff (c,ipp);
00375 Mc->give_area (eid,a);
00376 Mc->give_mominer (eid,ixyz);
00377 Mc->give_shearcoeff (eid,beta);
00378 Mc->give_vectorlcs (eid, vec);
00379 beam_transf_matrix (tran,dl, vec,x,y,z,eid);
00380 e=Mm->eliso[imat].e; g=e/(1.+2.*Mm->eliso[imat].nu);
00381
00382
00383 c1[0]=c[0][0];
00384 c1[1]=c[1][1];
00385 c1[2]=c[2][2];
00386 c2[0]=c[3][3];
00387 c2[1]=c[4][4];
00388 c2[2]=c[5][5];
00389
00390 k=0;
00391 for (i=0;i<intordsm[0][0];i++){
00392 Mm->computenlstresses (ipp);
00393 Mm->givestress (lcid,ipp,eps);
00394 for (j=0;j<6;j++){ rl[k]=eps[j]; k++;}
00395 ipp++;
00396 }
00397
00398 ll=dl*dl;
00399 f[0]= c1[0]*dl/3*rl[0]+ c1[0]*dl/6*rl[6];
00400 f[6]= c1[0]*dl/6*rl[0]+ c1[0]*dl/3*rl[6];
00401 g2=(c2[1]+c2[2])/2.;
00402 a2=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt((c1[1]+c1[2])/2.*g2) );
00403 a4=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt(g2*g2*g2/((c1[1]+c1[2])/2.)) );
00404 f[3]= (a2/3.*dl+a4/dl)*rl[3]+(a2/6.*dl-a4/dl)*rl(9);
00405 f[9]= (a2/6.*dl-a4/dl)*rl[3]+(a2/3.*dl+a4/dl)*rl(9);
00406
00407
00408 if(beta[0]<Mp->zero) gy=0.0;
00409 else gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
00410
00411 if(beta[1]<Mp->zero) gz=0.;
00412 else gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
00413
00414
00415 g2=gy*gy;
00416 a2=bPod*(c1[2]+sqrt(c1[2]*c2[2])*2./bPod)*dl/(1.+2.*gy)/(1.+2.*gy);
00417 a4=bPod*(c2[2]+sqrt(c2[2]*c2[2]*c2[2]/c1[2])/bPod);
00418
00419 f[2] = (a2*(4.*g2/3.+1.4*gy+13./35.)+a4*4*(g2+gy+0.3))*rl[2]
00420 +(a2*(2./3.*g2+0.6*gy+9./70.)-a4*4*(g2+gy+0.3))*rl[8]
00421 -(a2*(g2/6.+11.*gy/60.+11./210.)+a4/10.)*dl*rl[4]
00422 +(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl*rl[10];
00423
00424
00425
00426 f[4] = -(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl*rl[2]
00427 +(a2*(g2/30.+gy/30.+1./105.)+a4*(g2+gy+0.4)/3.)*ll*rl[4]
00428 -(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl*rl[8]
00429 -(a2*(g2/30.+gy/30.+1./140.)+a4*(g2+gy+0.1)/3.)*ll*rl[10];
00430 f[8] = (a2*(2./3.*g2+0.6*gy+9./70.)-a4*4*(g2+gy+0.3))*rl[2]
00431 +(a2*(4.*g2/3.+1.4*gy+13./35.)+a4*4*(g2+gy+0.3))*rl[8]
00432 -(a2*(g2+0.9*gy+13./70.)/6.-a4/10.)*dl*rl[4]
00433 +(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl*rl[10];
00434
00435
00436
00437 f[10]= +(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl*rl[2]
00438 -(a2*(g2/30.+gy/30.+1./140.)+a4*(g2+gy+0.1)/3.)*ll*rl[4]
00439 +(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl*rl[8]
00440 +(a2*(g2/30.+gy/30.+1./105.)+a4*(g2+gy+0.4)/3.)*ll*rl[10];
00441
00442
00443 g2=gz*gz;
00444 a2=bPod*(c1[1]+sqrt(c1[1]*c2[1])*2./bPod)*dl/(1.+2.*gz)/(1.+2.*gz);
00445 a4=bPod*(c2[1]+sqrt(c2[1]*c2[1]*c2[1]/c1[1])/bPod);
00446
00447
00448
00449 f[1] = (a2*(4.*g2/3.+1.4*gz+13./35.)+a4*4*(g2+gz+0.3))*rl[1]
00450 +(a2*(2./3.*g2+0.6*gz+9./70.)-a4*4*(g2+gz+0.3))*rl[7]
00451 -(a2*(g2/6.+11.*gz/60.+11./210.)+a4/10.)*dl*rl[5]
00452 +(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl*rl[11];
00453
00454
00455 f[5] = -(a2*(g2/6.+11.*gz/60.+11./210.)+ a4/10. )*dl*rl[1]
00456 +(a2*(g2/30.+gz/30.+1./105.)+a4*(g2+gz+0.4)/3.)*ll*rl[5]
00457 -(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl*rl[7]
00458 -(a2*(g2/30.+gz/30.+1./140.)+a4*(g2+gz+0.1)/3.)*ll*rl[11];
00459
00460
00461 f[7] = (a2*(2./3.*g2+0.6*gz+9./70.)-a4*4*(g2+gz+0.3))*rl[1]
00462 +(a2*(4.*g2/3.+1.4*gz+13./35.)+a4*4*(g2+gz+0.3))*rl[7]
00463 -(a2*(g2+0.9*gz+13./70.)/6.-a4/10.)*dl*rl[5]
00464 +(a2*(g2/6.+11.*gz/60.+11./210.)+ a4/10. )*dl*rl[11];
00465
00466
00467
00468
00469 f[11]= +(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl*rl[1]
00470 -(a2*(g2/30.+gz/30.+1./140.)+a4*(g2+gz+0.1)/3.)*ll*rl[5]
00471 +(a2*(g2/6.+11.*gz/60.+11./210.)+ a4/10. )*dl*rl[7]
00472 +(a2*(g2/30.+gz/30.+1./105.)+a4*(g2+gz+0.4)/3.)*ll*rl[11];
00473
00474
00475 for (k=0;k<ndofe;k++){
00476 ifor[k]+=f[k];
00477 }
00478
00479
00480
00481
00482
00483
00484
00485 }
00486 void soilbeam::internal_forces1 (long lcid,long eid,long ri,long ci,vector &ifor)
00487 {
00488 long i,j,k,integr=2,ipp;
00489 double dl,ll,a,ixyz[3],beta[2],g2,gy,gz,a2x,a2y,a2z,e,g,s;
00490 ivector cn(ndofe),nodes(nne);
00491 vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),f(ndofe),fx(tncomp),w(integr),gp(integr),eps(6);
00492 matrix n(tncomp,ndofe),c(tncomp,tncomp),d(tncomp,tncomp),tran(3,3),r(tncomp,2);
00493
00494 Mt->give_elemnodes (eid,nodes);
00495 Mt->give_node_coord3d (x,y,z,eid);
00496
00497 ipp=Mt->elements[eid].ipp[ri][ci];
00498 Mm->matstiff (c,ipp);
00499 Mm->elmatstiff (d,ipp);
00500 Mc->give_area (eid,a);
00501 Mc->give_mominer (eid,ixyz);
00502 Mc->give_shearcoeff (eid,beta);
00503 Mc->give_vectorlcs (eid, vec);
00504 beam_transf_matrix (tran,dl, vec,x,y,z,eid);
00505 e=d[0][0]; g=d[1][1]; ll=dl*dl;
00506 c1[0]=c[0][0];
00507 c1[1]=c[1][1];
00508 c1[2]=c[2][2];
00509 c2[0]=c[3][3];
00510 c2[1]=c[4][4];
00511 c2[2]=c[5][5];
00512
00513 k=0;
00514 for (i=0;i<intordsm[0][0];i++){
00515 Mm->computenlstresses (ipp);
00516 Mm->givestress (lcid,ipp,eps);
00517 for (j=0;j<6;j++){ rl[k]=eps[j]; k++; }
00518 ipp++;
00519 }
00520
00521
00522 if(beta[0]<Mp->zero) gy=0.0;
00523 else gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
00524
00525 if(beta[1]<Mp->zero) gz=0.;
00526 else gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
00527 g2=(c2[1]+c2[2])/2.;
00528 a2x=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt((c1[1]+c1[2])/2.*g2) );
00529 g2=gy*gy;
00530 a2y=bPod*(c1[2]+sqrt(c1[2]*c2[2])*2./bPod)*dl/(1.+2.*gy)/(1.+2.*gy);
00531 g2=gz*gz;
00532 a2z=bPod*(c1[1]+sqrt(c1[1]*c2[1])*2./bPod)*dl/(1.+2.*gz)/(1.+2.*gz);
00533
00534 s=0.;
00535 for (i=0;i<2;i++){
00536 geom_matrix (n,s, dl, gy, gz);
00537 fx[0]=( n[0][0]*rl[0]+n[0][6]*rl[6] )*c1[0];
00538 fx[3]=( n[3][3]*rl[3]+n[3][9]*rl[9] )*a2x;
00539 fx[1]=( n[1][1]*rl[1]+n[1][5]*rl[5]+n[1][7]*rl[7]+n[1][11]*rl[11] )*a2y;
00540 fx[5]=( n[5][1]*rl[1]+n[5][5]*rl[5]+n[5][7]*rl[7]+n[5][11]*rl[11] )*a2y;
00541 fx[2]=( n[2][2]*rl[2]+n[2][4]*rl[4]+n[2][8]*rl[8]+n[2][10]*rl[10] )*a2z;
00542 fx[4]=( n[4][2]*rl[2]+n[4][4]*rl[4]+n[4][8]*rl[8]+n[4][10]*rl[10] )*a2z;
00543 s=dl;
00544 }
00545
00546
00547
00548
00549
00550 }