00001 #include "beamgen3d.h"
00002 #include "global.h"
00003 #include "globmat.h"
00004 #include "genfile.h"
00005 #include "element.h"
00006 #include "node.h"
00007 #include "intpoints.h"
00008 #include <math.h>
00009
00010
00011 beamgen3d::beamgen3d (void)
00012 {
00013 long i,j;
00014
00015 nne=2; ndofe=14; tncomp=8; napfun=8;
00016 nb=1; intordmm=4; intordism=2; ssst=spacebeam;
00017
00018 ncomp = new long [nb];
00019 ncomp[0]=1;
00020
00021 cncomp = new long [nb];
00022 cncomp[0]=0;
00023
00024 nip = new long* [nb];
00025 intordsm = new long* [nb];
00026 for (i=0;i<nb;i++){
00027 nip[i] = new long [nb];
00028 intordsm[i] = new long [nb];
00029 }
00030
00031
00032
00033 nip[0][0]=2;
00034 intordsm[0][0]=1;
00035
00036 tnip=0;
00037 for (i=0;i<nb;i++){
00038 for (j=0;j<nb;j++){
00039 tnip+=nip[i][j];
00040 }
00041 }
00042 }
00043
00044 beamgen3d::~beamgen3d (void)
00045 {
00046 long i;
00047
00048 for (i=0;i<nb;i++){
00049 delete [] nip[i];
00050 }
00051 delete [] nip;
00052 }
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 void beamgen3d::transf_matrix (ivector &nodes,matrix &tmat)
00064 {
00065 long i,i6,n,m;
00066
00067 fillm (0.0,tmat);
00068
00069 n=nodes.n;
00070 m=tmat.m;
00071 for (i=0;i<m;i++){
00072 tmat[i][i]=1.0;
00073 }
00074
00075 for (i=0;i<n;i++){
00076 if (Mt->nodes[nodes[i]].transf>0){
00077 i6=i*7;
00078 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];
00079 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];
00080 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];
00081 i6=i*7+3;
00082 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];
00083 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];
00084 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];
00085 tmat[i6+3][i6] = 1.;
00086 }
00087 }
00088 }
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 void beamgen3d::beam_transf_matrix (matrix &tmat,double &dl,vector &vec,vector &x,vector &y,vector &z,long eid)
00106 {
00107 double c;
00108
00109 fillm (0.0,tmat);
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
00116
00117 dl=sqrt((tmat[0][0]*tmat[0][0]+tmat[1][0]*tmat[1][0]+tmat[2][0]*tmat[2][0]));
00118
00119 if (dl<Mp->zero){
00120 fprintf (stderr,"\n\n zero length of the %ld beamgen3d element",eid);
00121 fprintf (stderr,"\n in function beamgen3d::beam_transf_matrix (file %s, line %d).\n",__FILE__,__LINE__);
00122 }
00123
00124
00125 tmat[0][0]=tmat[0][0]/dl;
00126 tmat[1][0]=tmat[1][0]/dl;
00127 tmat[2][0]=tmat[2][0]/dl;
00128
00129
00130 tmat[0][1]=vec[1]*tmat[2][0]-vec[2]*tmat[1][0];
00131 tmat[1][1]=vec[2]*tmat[0][0]-vec[0]*tmat[2][0];
00132 tmat[2][1]=vec[0]*tmat[1][0]-vec[1]*tmat[0][0];
00133
00134
00135 c=sqrt((tmat[0][1]*tmat[0][1]+tmat[1][1]*tmat[1][1]+tmat[2][1]*tmat[2][1]));
00136
00137 if (c<Mp->zero){
00138 vec[0]=0.0;
00139 vec[1]=1.0;
00140 vec[2]=0.0;
00141
00142
00143 tmat[0][2]=tmat[1][0]*vec[2]-tmat[2][0]*vec[1];
00144 tmat[1][2]=tmat[2][0]*vec[0]-tmat[0][0]*vec[2];
00145 tmat[2][2]=tmat[0][0]*vec[1]-tmat[1][0]*vec[0];
00146
00147
00148 c=sqrt((tmat[0][2]*tmat[0][2]+tmat[1][2]*tmat[1][2]+tmat[2][2]*tmat[2][2]));
00149
00150
00151 tmat[0][2]=tmat[0][2]/c;
00152 tmat[1][2]=tmat[1][2]/c;
00153 tmat[2][2]=tmat[2][2]/c;
00154
00155 if (c<Mp->zero){
00156 fprintf (stderr,"\n\n zero z base vector of the %ld beamgen3d element",eid);
00157 fprintf (stderr,"\n in function beamel2d::beam_transf_matrix.\n");
00158 }
00159
00160
00161 tmat[0][1]=tmat[1][2]*tmat[2][0]-tmat[2][2]*tmat[1][0];
00162 tmat[1][1]=tmat[2][2]*tmat[0][0]-tmat[0][2]*tmat[2][0];
00163 tmat[2][1]=tmat[0][2]*tmat[1][0]-tmat[1][2]*tmat[0][0];
00164
00165
00166 c=sqrt((tmat[0][1]*tmat[0][1]+tmat[1][1]*tmat[1][1]+tmat[2][1]*tmat[2][1]));
00167
00168 if (c<Mp->zero){
00169 fprintf (stderr,"\n\n zero y base vector of the %ld beamgen3d element",eid);
00170 fprintf (stderr,"\n in function beamgen3d::beam_transf_matrix (file %s, line %d).\n",__FILE__,__LINE__);
00171 }
00172
00173
00174 tmat[0][1]=tmat[0][1]/c;
00175 tmat[1][1]=tmat[1][1]/c;
00176 tmat[2][1]=tmat[2][1]/c;
00177
00178 }
00179 else{
00180
00181
00182 tmat[0][1]=tmat[0][1]/c;
00183 tmat[1][1]=tmat[1][1]/c;
00184 tmat[2][1]=tmat[2][1]/c;
00185
00186
00187 tmat[0][2]=tmat[1][0]*tmat[2][1]-tmat[2][0]*tmat[1][1];
00188 tmat[1][2]=tmat[2][0]*tmat[0][1]-tmat[0][0]*tmat[2][1];
00189 tmat[2][2]=tmat[0][0]*tmat[1][1]-tmat[1][0]*tmat[0][1];
00190
00191
00192 c=sqrt((tmat[0][2]*tmat[0][2]+tmat[1][2]*tmat[1][2]+tmat[2][2]*tmat[2][2]));
00193
00194 if (c<Mp->zero){
00195 fprintf (stderr,"\n\n zero z base vector of the %ld beamgen3d element",eid);
00196 fprintf (stderr,"\n in function beamel2d::beam_transf_matrix.\n");
00197 }
00198
00199
00200 tmat[0][2]=tmat[0][2]/c;
00201 tmat[1][2]=tmat[1][2]/c;
00202 tmat[2][2]=tmat[2][2]/c;
00203 }
00204 }
00205 void beamgen3d::ck_matrix (matrix &ck, double s,double c,double eh,double dl)
00206 {
00207 double eh1,c1,cjm;
00208
00209 fillm (0.0,ck);
00210 eh1=eh*dl;
00211
00212
00213 c1=c-1.;
00214 cjm=c1*c1-(s-eh1)*s;
00215
00216 ck[0][3]=s/cjm;
00217 ck[0][2]=-ck[0][3]*c1/s;
00218 ck[0][1]=-eh*ck[0][3];
00219 ck[0][0]=1.-ck[0][2];
00220
00221 ck[1][3]= (eh1*s-c1)/eh/cjm;
00222 ck[1][2]= -(1./eh+ck[1][3]*c1)/s;
00223 ck[1][1]= 1.-eh*ck[1][3];
00224 ck[1][0]=-ck[1][2];
00225
00226 ck[2][3]=-ck[0][3];
00227 ck[2][2]=-ck[2][3]*c1/s;
00228 ck[2][1]=-eh*ck[2][3];
00229 ck[2][0]=-ck[2][2];
00230
00231 ck[3][3]= c1/eh/cjm;
00232 ck[3][2]=-ck[3][3]*(s-eh1)/c1;
00233 ck[3][1]=-eh*ck[3][3];
00234 ck[3][0]=-ck[3][2];
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 void beamgen3d::geom_matrix (matrix &n,double s,double dl,double gy,double gz)
00260 {
00261 double aj1,ll;
00262
00263 fillm (0.0,n);
00264 ll=dl*dl;
00265
00266
00267 n[0][0]=-1./dl;
00268 n[0][6]= 1./dl;
00269
00270 aj1=1./(1.+2*gy);
00271 n[4][2] = aj1*(6.-12.*s)/ll;
00272 n[4][4] = aj1*(-2.*(2.+gy)+6.*s)/dl;
00273 n[4][8] =-n[4][2];
00274 n[4][10]= aj1*(-2.*(1.-gy)+6.*s)/dl;
00275
00276 n[2][2] =-2.*gy/dl/(1.+2.*gy);
00277 n[2][4] = gy/(1.+2.*gy);
00278 n[2][8] =-n[2][2];
00279 n[2][10]= n[2][4];
00280
00281 aj1=1./(1.+2.*gz);
00282 n[5][1] =-aj1*( 6.-12.*s)/ll;
00283 n[5][5] = aj1*(-2.*(2.+gz)+6.*s)/dl;
00284 n[5][7] =-n[5][1];
00285 n[5][11]= aj1*(-2.*(1.-gz)+6.*s)/dl;
00286
00287 n[1][1] =-2.*gz/dl/(1.+2.*gz);
00288 n[1][5] =-gz/(1.+2.*gz);
00289 n[1][7] =-n[1][1];
00290 n[1][11]= n[1][5];
00291
00292 n[3][3] =-1./dl;
00293 n[3][9] = 1./dl;
00294 }
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311 void beamgen3d::bf_matrix (matrix &n,double s,double dl,double gy,double gz)
00312 {
00313 double ll,aj1;
00314 ll=dl*dl;
00315
00316 n[0][0] = 1.-s;
00317 n[0][6] = s;
00318
00319 aj1=1./(1.+2.*gy);
00320 n[2][2] =aj1*(1.+2.*gy-2.*gy*s-3.*s*s+2*s*s*s);
00321 n[2][4] =aj1*(-(1.+gy)*s+(2.+gy)*s*s-s*s*s)*dl;
00322 n[2][8] =aj1*(2.*gy*s+3.*s*s-2.*s*s*s);
00323 n[2][10]=aj1*(gy*s+(1.-gy)*s*s-s*s*s)*dl;
00324
00325 n[4][2] = aj1*(6.*s-6.*s*s)/dl;
00326 n[4][4] = aj1*(1.+2.*gy-2.*(2.+gy)*s+3.*s*s);
00327 n[4][8] =-n[4][2];
00328 n[4][10]= aj1*(-2.*(1.-gy)*s+3.*s*s);
00329
00330
00331 aj1=1./(1.+2.*gz);
00332 n[1][1] = aj1*(1.+2.*gz-2.*gz*s-3.*s*s+2.*s*s*s);
00333 n[1][5] =-aj1*(-(1.+gz)*s+(2.+gz)*s*s-s*s*s)*dl;
00334 n[1][7] = aj1*(2.*gz*s+3.*s*s-2.*s*s*s);
00335 n[1][11]=-aj1*(gz*s+(1.-gz)*s*s-s*s*s)*dl;
00336
00337 n[5][1] =-aj1*(6.*s-6.*s*s)/dl;
00338 n[5][5] = aj1*(1.+2.*gz-2.*(2.+gz)*s+3.*s*s);
00339 n[5][7] =-n[5][1];
00340 n[5][11]= aj1*(-2.*(1.-gz)*s+3.*s*s);
00341
00342 n[3][3] = 1.-s;
00343 n[3][9] = s;
00344
00345 }
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360 void beamgen3d::res_stiffness_matrix (long eid,matrix &sm)
00361 {
00362 stiffness_matrix (eid,0,0,sm);
00363
00364 long i;
00365 for (i=0;i<sm.m;i++){
00366 if (sm[i][i]<Mp->zero){
00367 fprintf (stderr,"\n\n nulovy diagonalni prvek v matici tuhosti jednoho prvku (file %s, line %d)",__FILE__,__LINE__);
00368 }
00369 }
00370 }
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 void beamgen3d::stiffness_matrix (long eid,long ri,long ci,matrix &sm)
00382 {
00383 long ipp,i,transf;
00384 double e,g,a,*ixyz,*ioyz,*beta,dl,ll,gy,gz,aj1, integn4;
00385 ivector nodes(nne);
00386 vector vec(3),x(nne),y(nne),z(nne),integn1(4),integn2(3),integn3(2);
00387 matrix d(tncomp,tncomp),tran(3,3);
00388
00389 ixyz = new double [3];
00390 ioyz = new double [3];
00391 beta = new double [2];
00392
00393 Mt->give_elemnodes (eid,nodes);
00394 Mt->give_node_coord3d (x,y,z,eid);
00395 Mc->give_vectorlcs (eid, vec);
00396 beam_transf_matrix (tran,dl, vec,x,y,z,eid);
00397 ll=dl*dl;
00398
00399 fillm (0.0,sm);
00400 ipp=Mt->elements[eid].ipp[ri][ci];
00401 Mm->matstiff (d,ipp);
00402 Mc->give_area (eid,a);
00403 Mc->give_mominer (eid,ixyz);
00404 Mc->give_shearcoeff (eid,beta);
00405
00406
00407 ioyz[0]=4.2692e-6 ; ioyz[1]=-2.6687e-5; ioyz[2]=0.0; ixyz[0]=2.646e-7; ixyz[1]=2.135e-4; ixyz[2]=3.333e-5;
00408 ioyz[0]=1.1924e-6 ; ioyz[1]=0.0; ioyz[2]=7.5625e-6; ixyz[0]=2.333e-7; ixyz[1]=3.048e-5; ixyz[2]=6.25e-5;
00409
00410
00411
00412
00413 e=d[0][0]; g=d[1][1];
00414
00415
00416
00417 sm[0][0]= e*a/dl; sm[0][7]= -sm[0][0]; sm[7][7]= sm[0][0];
00418
00419
00420 if(beta[0]<Mp->zero) gy=0.0;
00421 else gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
00422 if(beta[1]<Mp->zero) gz=0.;
00423 else gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
00424 gy=0.;gz=0.;
00425 aj1=ixyz[1]*e/(1.+2*gy);
00426 integn1[0]=12./ll/dl*aj1; integn1[1]= 6./ll*aj1; integn1[2]=-12./ll/dl*aj1; integn1[3]= 6./ll*aj1;
00427 integn2[0]=2.*(2.+gy)/dl*aj1; integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-gy)/dl*aj1;
00428 integn3[0]=12./ll/dl*aj1; integn3[1]=-6./ll*aj1;
00429 integn4 =2.*(2.+gy)/dl*aj1;
00430 sm[2][2] = integn1[0];
00431 sm[2][4] = -integn1[1];
00432 sm[2][9] = integn1[2];
00433 sm[2][11]= -integn1[3];
00434 sm[4][4] = integn2[0];
00435 sm[4][9] = integn2[1];
00436 sm[4][11]= -integn2[2];
00437 sm[9][9] = integn3[0];
00438 sm[9][11]= -integn3[1];
00439 sm[11][11]= integn4;
00440
00441 aj1=ixyz[2]*e/(1.+2.*gz);
00442 integn1[0]=12./ll/dl*aj1; integn1[1]= 6./ll*aj1; integn1[2]=-12./ll/dl*aj1; integn1[3]=6./ll*aj1;
00443 integn2[0]=2.*(2.+gz)/dl*aj1; integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-gz)/dl*aj1;
00444 integn3[0]=12./ll/dl*aj1; integn3[1]=-6./ll*aj1;
00445 integn4 =2.*(2.+gz)/dl*aj1;
00446 sm[1][1] = integn1[0];
00447 sm[1][5] = integn1[1];
00448 sm[1][8] = integn1[2];
00449 sm[1][12]= integn1[3];
00450 sm[5][5] = integn2[0];
00451 sm[5][8] = integn2[1];
00452 sm[5][12]= integn2[2];
00453 sm[8][8] = integn3[0];
00454 sm[8][12]= integn3[1];
00455 sm[12][12]= integn4;
00456
00457 stiffness_matrixtor (sm, dl,e,g,gy,gz,ixyz,ioyz);
00458
00459 for (i=0;i<ndofe;i++){for (int j=i+1;j<ndofe;j++){sm[j][i]=sm[i][j];}}
00460
00461
00462
00463
00464
00465 transf = Mt->locsystems (nodes);
00466 if (transf>0){
00467 matrix tmat (ndofe,ndofe);
00468 transf_matrix (nodes,tmat);
00469 glmatrixtransf (sm,tmat);
00470 }
00471
00472 fprintf (Out,"\n\n MT prvek beam cislo %ld",eid);
00473 for (i=0;i<ndofe;i++){
00474 fprintf (Out,"\n %4ld",i);
00475 for (int j=0;j<ndofe;j++){
00476 fprintf (Out," %15.5e",sm[i][j]);
00477 }
00478 }
00479
00480 delete [] ixyz;
00481 delete [] ioyz;
00482 delete [] beta;
00483 }
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 void beamgen3d::stiffness_matrixtor (matrix &sm, double dl,double e,double g,double gy,double gz,double *ixyz,double *ioyz)
00499 {
00500 int i,i1,i2,j;
00501 double eh,eh1,s,c,c1,s1,s2,c2,cs,pom1,pom2,pom3,xc,xs;
00502 vector sioyz(2),sm1(16),sm2(16);
00503 matrix ck(4,4);
00504
00505
00506
00507
00508
00509 eh=pow(g*ixyz[0]/e/ioyz[0],0.5);
00510
00511 eh1=eh*dl;
00512 s=sinh(eh1); c=cosh(eh1);
00513 ck_matrix (ck, s,c,eh,dl);
00514
00515
00516 c1=s/eh; s1=(c-1.)/eh; xc = (dl*s-s1)/eh; xs = (dl*c-c1)/eh; c2=s*c/2/eh+dl/2; cs=s*s/2/eh; s2=s*c/2/eh-dl/2;
00517 pom1=e*ioyz[0]*eh*eh*eh*eh;
00518 pom2=g*ixyz[0];
00519
00520 pom3=pom1/eh/eh;
00521 for (i=0;i<4;i++){
00522 i1=4*i-i*(i+1)/2;
00523 for (j=i;j<4;j++){
00524 i2=j+i1;
00525 sm1[i2]=(ck[i][2]*ck[j][2]*c2+ck[i][3]*ck[j][3]*s2+
00526 (ck[i][3]*ck[j][2]+ck[i][2]*ck[j][3])*cs)*pom1;
00527 sm1[i2]=sm1[i2]+(ck[i][1]*ck[j][1]*dl+eh*(ck[i][1]*ck[j][2]+
00528 ck[i][2]*ck[j][1])*s1+eh*(ck[i][1]*ck[j][3]+ck[i][3]*ck[j][1])*c1)*pom2;
00529 sm1[i2]=sm1[i2]+((ck[i][2]*ck[j][3]+ck[i][3]*ck[j][2])*cs+
00530 ck[i][2]*ck[j][2]*s2+ck[i][3]*ck[j][3]*c2)*eh*eh*pom2;
00531 }
00532 }
00533
00534 sm[3][3] = sm1[0] ;
00535 sm[3][6] = sm1[1] ;
00536 sm[3][10] = sm1[2] ;
00537 sm[3][13] = sm1[3] ;
00538 sm[6][6] = sm1[4] ;
00539 sm[6][10] = sm1[5] ;
00540 sm[6][13] = sm1[6] ;
00541 sm[10][10]= sm1[7] ;
00542 sm[10][13]= sm1[8] ;
00543 sm[13][13]= sm1[9] ;
00544
00545
00546
00547
00548
00549 for (i=0;i<4;i++){
00550 sm1[i] = 6.*eh*eh/dl/dl*(-(ck[i][2]*c1+ck[i][3]*s1)+
00551 2./dl*(ck[i][2]*xc+ck[i][3]*xs) );
00552 sm1[i+4] =-2.*eh*eh/dl*(-(2.+gy)*(ck[i][2]*c1+ck[i][3]*s1)+
00553 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
00554 sm2[i+4] =-2.*eh*eh/dl*(-(2.+gz)*(ck[i][2]*c1+ck[i][3]*s1)+
00555 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
00556 sm1[i+8] =-sm1[i];
00557 sm1[i+12]=-2.*eh*eh/dl*(-(1.-gy)*(ck[i][2]*c1+ck[i][3]*s1)+
00558 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
00559 sm2[i+12]=-2.*eh*eh/dl*(-(1.-gz)*(ck[i][2]*c1+ck[i][3]*s1)+
00560 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
00561 }
00562
00563 sioyz[0]=ioyz[1]*e/(1.+2.*gy);
00564 sioyz[1]=ioyz[2]*e/(1.+2.*gz);
00565 for (i=0;i<2;i++){
00566 i1=1-i;
00567 i2=i+1;
00568 sm[i2][3] = sm1[0]*sioyz[i1];
00569 sm[i2][6] = sm1[1]*sioyz[i1];
00570 sm[i2][10]= sm1[2]*sioyz[i1];
00571 sm[i2][13]= sm1[3]*sioyz[i1];
00572 sm[3][4] = sm1[4]*sioyz[0];
00573 sm[3][5] =-sm2[4]*sioyz[1];
00574 sm[3][i+8]= sm1[8]*sioyz[i1];
00575 sm[3][11] = sm1[12]*sioyz[0];
00576 sm[3][12] =-sm2[12]*sioyz[1];
00577
00578 sm[4][6] = sm1[5]*sioyz[0];
00579 sm[5][6] =-sm2[5]*sioyz[1];
00580 sm[4][10] = sm1[6]*sioyz[0];
00581 sm[5][10] =-sm2[6]*sioyz[1];
00582 sm[4][13] = sm1[7]*sioyz[0];
00583 sm[5][13] =-sm2[7]*sioyz[1];
00584 sm[6][i+8]= sm1[9]*sioyz[i1];
00585 sm[6][11] = sm1[13]*sioyz[0];
00586 sm[6][12] =-sm2[13]*sioyz[1];
00587 i2=i+8;
00588 sm[i2][10]= sm1[10]*sioyz[i1];
00589 sm[i2][13]= sm1[11]*sioyz[i1];
00590 sm[10][11]= sm1[14]*sioyz[0];
00591 sm[10][12]=-sm2[14]*sioyz[1];
00592
00593 sm[11][13]= sm1[15]*sioyz[0];
00594 sm[12][13]=-sm2[15]*sioyz[1];
00595 }
00596
00597 }
00598
00599
00600
00601
00602 void beamgen3d::stiffness_matrixtor1 (matrix &sm, double dl,double e,double g,double gy,double gz,double *ixyz,double *ioyz)
00603 {
00604 double aj1,ll,integn4;
00605 vector sioyz(2),integn1(4),integn2(3),integn3(2);
00606
00607
00608
00609
00610
00611 ll=dl*dl; aj1=e*ioyz[0];
00612 integn1[0]=12./ll/dl*aj1; integn1[1]= 6./ll*aj1; integn1[2]=-12./ll/dl*aj1; integn1[3]= 6./ll*aj1;
00613 integn2[0]=2.*(2.+0.)/dl*aj1; integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-0.)/dl*aj1;
00614 integn3[0]=12./ll/dl*aj1; integn3[1]=-6./ll*aj1;
00615 integn4=2.*(2.+0.)/dl*aj1;
00616
00617 sm[3][3] = integn1[0]+g*1.2*ixyz[0]/dl;
00618 sm[3][6] = integn1[1]+g*0.1*ixyz[0];
00619 sm[3][10] = integn1[2]-g*1.2*ixyz[0]/dl;
00620 sm[3][13] = integn1[3]+g*0.1*ixyz[0];
00621 sm[6][6] = integn2[0]+g*4./30.*ixyz[0]*dl;
00622 sm[6][10] = integn2[1]-g*0.1*ixyz[0];
00623 sm[6][13] = integn2[2]-g*1./30.*ixyz[0]*dl;
00624 sm[10][10]= integn3[0]+g*1.2*ixyz[0]/dl;
00625 sm[10][13]= integn3[1]-g*0.1*ixyz[0];
00626 sm[13][13]= integn4 +g*4./30.*ixyz[0]*dl;
00627
00628
00629 aj1=ioyz[1]*e/(1.+2.*gy);
00630 integn1[0]=12./ll/dl*aj1; integn1[1]= 6./ll*aj1; integn1[2]=-12./ll/dl*aj1; integn1[3]= 6./ll*aj1;
00631 integn2[0]=2.*(2.+gy)/dl*aj1; integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-gy)/dl*aj1;
00632 integn3[0]=12./ll/dl*aj1; integn3[1]=-6./ll*aj1;
00633 integn4=2.*(2.+gy)/dl*aj1;
00634
00635 sm[2][3] = integn1[0];
00636 sm[2][6] = integn1[1];
00637 sm[2][10]= integn1[2];
00638 sm[2][13]= integn1[3];
00639 sm[3][4] = -integn1[1];
00640 sm[3][8]=integn1[2];
00641 sm[3][11] = -integn1[3];
00642 sm[4][6] = -integn2[0];
00643 sm[4][10] = -integn2[1];
00644 sm[4][13] = -integn2[2];
00645 sm[6][9]= integn2[1];
00646 sm[6][11] = -integn2[2];
00647 sm[9][10]= integn3[0];
00648 sm[9][13]= integn3[1];
00649 sm[10][11]= -integn3[1];
00650 sm[11][13]= -integn4;
00651
00652
00653 aj1=ioyz[2]*e/(1.+2.*gz);
00654 integn1[0]=12./ll/dl*aj1; integn1[1]= 6./ll*aj1; integn1[2]=-12./ll/dl*aj1; integn1[3]= 6./ll*aj1;
00655 integn2[0]=2.*(2.+gz)/dl*aj1; integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-gz)/dl*aj1;
00656 integn3[0]=12./ll/dl*aj1; integn3[1]=-6./ll*aj1;
00657 integn4=2.*(2.+gz)/dl*aj1;
00658 sm[1][3] = integn1[0];
00659 sm[1][6] = integn1[1];
00660 sm[1][10]= integn1[2];
00661 sm[1][13]= integn1[3];
00662 sm[3][5] = integn1[1];
00663 sm[3][8]=integn1[2];
00664 sm[3][12] = integn1[3];
00665 sm[5][6] = integn2[0];
00666 sm[5][10] = integn2[1];
00667 sm[5][13] = integn2[2];
00668 sm[6][8]= integn2[1];
00669 sm[6][12] = integn2[2];
00670 sm[8][10]= integn3[0];
00671 sm[8][13]= integn3[1];
00672 sm[10][12]= integn2[1];
00673 sm[12][13]= integn4;
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706 }
00707
00708
00709
00710
00711 void beamgen3d::stiffness_matrixtor2 (matrix &sm, double dl,double e,double g,double gy,double gz,double *ixyz,double *ioyz,double *iro)
00712 {
00713 double kakr,ajy,ajz,ajk;
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724 kakr =6.*ioyz[0]*e/iro[0]/g/dl/dl;
00725
00726 ajy= (1.+2.*gy);
00727 ajk= (1.+2.*kakr);
00728 ajy=1.;
00729 sm[2][3] = (12.*ioyz[1]*e/dl/dl/dl+4.*kakr*gy*iro[1]*g/dl)/ajk/ajy;
00730 sm[3][3] = (12.*ioyz[0]*e/dl/dl/dl+(4.*kakr*kakr*iro[0]
00731 +ixyz[0]*(4.*kakr*(kakr+1.)+1.2) )*g/dl )/ajk/ajk;
00732 sm[3][4] = (-6.*ioyz[1]*e/dl/dl-2*kakr*gy*iro[1]*g)/ajk/ajy;
00733 sm[2][6] = sm[4][5];
00734 sm[3][6] = (-6.*ioyz[0]*e/dl/dl-(2.*kakr*kakr*iro[0]+iro[0]/10.)*g )/ajk/ajk;
00735 sm[4][6] = ((2.+kakr+gy+2.*kakr*gy)*2.*ioyz[1]*e/dl
00736 +kakr*gy*iro[1]*g*dl)/ajk/ajy;
00737 sm[6][6] = ((1.+(1.+kakr)*kakr)*4.*ioyz[0]*e/dl+(kakr*kakr*iro[0]
00738 +(kakr*(kakr+1)+0.4)/3.*ixyz[0])*g*dl)/ajk/ajk;
00739
00740 sm[3][9] = -sm[2][3]; sm[6][9] = -sm[3][4];
00741 sm[2][10] = -sm[2][3]; sm[3][10] = -sm[3][3];
00742 sm[4][10] = -sm[3][4]; sm[6][10] = -sm[3][6];
00743 sm[9][10] = sm[2][3]; sm[10][10]= sm[3][3];
00744
00745 sm[3][11] = sm[3][4];
00746 sm[6][11] = ((1.-kakr-gy-2.*kakr*gy)*2.*ioyz[1]*e/dl
00747 -kakr*gy*iro[1]*g*dl)/ajk/ajy;
00748 sm[10][11]= -sm[3][4]; sm[2][13] = sm[2][6];
00749 sm[3][13] = sm[3][6];
00750 sm[4][13] = ((1.-kakr-gy-2.*kakr*gy)*2.*ioyz[1]*e/dl
00751 +kakr*gy*iro[1]*g*dl)/ajk/ajy;
00752 sm[6][13] = ((2.-(4.+4.*kakr)*kakr)*ioyz[0]*e/dl+(kakr*kakr*iro[0]
00753 +(kakr*kakr-0.1)/3.*ixyz[0] )*g*dl)/ajk/ajk;
00754 sm[9][13] = -sm[2][6];sm[10][13] = -sm[3][6];
00755 sm[11][13] = sm[4][6];sm[13][13] = sm[6][6];
00756
00757
00758 ajz= (1.+2.*gz);
00759 ajz=1.;
00760 sm[1][3] = (12.*ioyz[2]*e/dl/dl/dl+4.*kakr*gz*iro[2]*g/dl)/ajk/ajz;
00761 sm[3][5] = (-6.*ioyz[2]*e/dl/dl-2.*kakr*gz*iro[2]*g)/ajk/ajz;
00762 sm[1][6] = sm[3][5];
00763 sm[5][6] = ((2.+kakr+gz+2.*kakr*gz)*2.*ioyz[2]*e/dl
00764 +kakr*gz*iro[2]*g*dl)/ajk/ajz;
00765
00766 sm[3][8] = -sm[1][3];sm[6][8] = -sm[3][5];
00767 sm[1][10] = -sm[1][3];sm[5][10] = -sm[3][5];
00768 sm[8][10] = sm[1][3];
00769
00770 sm[3][12] = sm[3][5];
00771 sm[6][12] = ((1.-kakr-gz-2.*kakr*gz)*2.*ioyz[2]*e/dl
00772 -kakr*gz*iro[2]*g*dl)/ajk/ajz;
00773 sm[10][12]= -sm[4][6];sm[2][14] = sm[2][7];
00774 sm[5][13] = ((1.-kakr-gz-2.*kakr*gz)*2.*ioyz[2]*e/dl
00775 +kakr*gz*iro[2]*g*dl)/ajk/ajz;
00776 sm[8][13] = -sm[1][6];
00777 sm[12][13]= sm[5][6];
00778
00779 }
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790 void beamgen3d::load_matrix (long eid,matrix &lm)
00791 {
00792 long i,ii,transf;
00793 double jac,a,ixyz[3],beta[2],e,g,gy,gz,l;
00794 ivector nodes (nne);
00795 vector vec(3),x(nne),y(nne),z(nne),w(intordmm),gp(intordmm),xl(2);
00796 matrix d(tncomp,tncomp),n(napfun,ndofe),tran(3,3);
00797
00798
00799 Mt->give_node_coord3d (x,y,z,eid);
00800
00801 Mc->give_vectorlcs (eid,vec);
00802
00803 beam_transf_matrix (tran,l,vec,x,y,z,eid);
00804
00805 Mc->give_area (eid,a);
00806
00807 Mc->give_shearcoeff (eid,beta);
00808
00809 Mc->give_mominer (eid,ixyz);
00810
00811 ii=Mt->elements[eid].ipp[0][0];
00812
00813 Mm->matstiff (d,ii);
00814
00815 e=d[0][0];
00816
00817 g=d[1][1];
00818
00819
00820 if(beta[0]<Mp->zero) gy=0.0;
00821 else gy=6.0*e*ixyz[1]/beta[0]/g/a/l;
00822
00823 if(beta[1]<Mp->zero) gz=0.;
00824 else gz=6.0*e*ixyz[2]/beta[1]/g/a/l;
00825
00826
00827
00828 xl[0]=0.0;
00829 xl[1]=l;
00830
00831
00832 Mt->give_elemnodes (eid,nodes);
00833
00834
00835 gauss_points (gp.a,w.a,intordmm);
00836
00837
00838 fillm (0.0,lm);
00839
00840 for (i=0;i<intordmm;i++){
00841
00842 jac_1d (jac,xl,gp[i]);
00843
00844
00845 bf_matrix (n,(gp[i]+1.0)/2.0,l,gy,gz);
00846
00847 jac*=a*w[i];
00848
00849
00850 nnj (lm.a,n.a,jac,n.m,n.n);
00851 }
00852
00853
00854
00855
00856
00857 transf = Mt->locsystems (nodes);
00858 if (transf>0){
00859 matrix tmat (ndofe,ndofe);
00860
00861 transf_matrix (nodes,tmat);
00862 glmatrixtransf (lm,tmat);
00863 }
00864 }
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877 void beamgen3d::nodal_displ (long eid,long lcid)
00878 {
00879 long ii,transf;
00880 double l;
00881 ivector nodes(nne);
00882 vector x(nne),y(nne),z(nne),vec(3),d(ndofe),f(ndofe);
00883 matrix tm(3,3);
00884
00885
00886 eldispl (lcid,eid,d.a);
00887
00888
00889 Mt->give_elemnodes (eid,nodes);
00890
00891
00892 transf = Mt->locsystems (nodes);
00893 if (transf>0){
00894 matrix tmat (ndofe,ndofe);
00895
00896 transf_matrix (nodes,tmat);
00897
00898 mtxv (tmat,d,f);
00899 copyv (f,d);
00900 }
00901
00902
00903 Mt->give_node_coord3d (x,y,z,eid);
00904
00905
00906 Mc->give_vectorlcs (eid,vec);
00907
00908
00909 beam_transf_matrix (tm,l,vec,x,y,z,eid);
00910
00911
00912
00913
00914
00915 ii=Mt->elements[eid].ipp[0][0];
00916
00917
00918 Mm->storestrain (lcid,ii,0,7,d);
00919 Mm->storestrain (lcid,ii+1,7,7,d);
00920
00921 }
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937 void beamgen3d::nodal_forces (long eid,long lcid)
00938 {
00939 long ii,j,transf;
00940 double l;
00941 ivector nodes(nne);
00942 vector x(nne),y(nne),z(nne),vec(3),d(ndofe),f(ndofe);
00943 matrix sm(ndofe,ndofe),tm(3,3);
00944
00945
00946 eldispl (lcid,eid,d.a);
00947
00948
00949 stiffness_matrix (eid,0,0,sm);
00950
00951
00952 mxv (sm,d,f);
00953 fprintf (Out,"\n\n MT prvek beam cislo %ld",lcid);
00954 fprintf (Out," \n u, v, w, fx, fy, fz, o \n ");
00955 for (j=0;j<ndofe;j++){
00956 fprintf (Out," %15.5e",d[j]);
00957 }
00958 fprintf (Out," \n N, Vy, Vz, Mx, My, Mz, B,\n ");
00959 for (j=0;j<ndofe;j++){
00960 fprintf (Out," %15.5e",f[j]);
00961 }
00962
00963
00964 Mt->give_elemnodes (eid,nodes);
00965
00966
00967 transf = Mt->locsystems (nodes);
00968 if (transf>0){
00969 matrix tmat (ndofe,ndofe);
00970
00971 transf_matrix (nodes,tmat);
00972
00973 mtxv (tmat,f,d);
00974 copyv (d,f);
00975 }
00976
00977
00978 Mt->give_node_coord3d (x,y,z,eid);
00979
00980
00981 Mc->give_vectorlcs (eid,vec);
00982
00983
00984 beam_transf_matrix (tm,l,vec,x,y,z,eid);
00985
00986 internal_forces (lcid, eid, 0, 0, f);
00987
00988
00989
00990
00991
00992 ii=Mt->elements[eid].ipp[0][0];
00993
00994
00995
00996 Mm->storestress (lcid,ii,0,6,f);
00997 Mm->storestress (lcid,ii+1,6,6,f);
00998
00999
01000 }
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012 void beamgen3d::res_internal_forces (long lcid,long eid,vector &ifor)
01013 {
01014 internal_forces1 (lcid, eid, 0, 0, ifor);
01015 }
01016
01017
01018 void beamgen3d::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
01019 {
01020 long i,j,integr=2,ipp,transf;
01021 double dl,ll,a,ixyz[3],ioyz[3],beta[2],gy,gz,kakr,e,g;
01022 double eh,ehx,ajy,ajz,ajk,xs, s,c;
01023 ivector nodes(nne);
01024 vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),f(ndofe),fx(tncomp),w(integr),gp(integr);
01025 vector dderno(4),dderny(4),ddernz(4);
01026 matrix n(tncomp,ndofe),d(tncomp,tncomp),tran(3,3),r(tncomp,2);
01027 matrix bb(tncomp,ndofe), ck(4,4);
01028
01029
01030
01031 Mt->give_elemnodes (eid,nodes);
01032 Mt->give_node_coord3d (x,y,z,eid);
01033 eldispl (lcid,eid,rl.a);
01034 Mc->give_vectorlcs (eid, vec);
01035 beam_transf_matrix (tran,dl, vec,x,y,z,eid);
01036 ll=dl*dl;
01037
01038
01039 fillm (0.0,r);
01040
01041
01042 transf = Mt->locsystems (nodes);
01043 if (transf>0){
01044 matrix tmat (ndofe,ndofe);
01045 transf_matrix (nodes,tmat);
01046
01047 lgvectortransf (rg,rl,tmat);
01048 }
01049 else{
01050 copyv (rl,rg);
01051 }
01052
01053
01054
01055
01056
01057
01058 ipp=Mt->elements[eid].ipp[0][0];
01059 Mm->matstiff (d,ipp);
01060 Mc->give_area (eid,a);
01061 Mc->give_mominer (eid,ixyz);
01062 Mc->give_shearcoeff (eid,beta);
01063
01064
01065 ioyz[0]=4.2692e-6 ; ioyz[1]=-2.6687e-5; ioyz[2]=0.0; ixyz[0]=2.646e-7; ixyz[1]=2.135e-4; ixyz[2]=3.333e-5;
01066 ioyz[0]=1.1924e-6 ; ioyz[1]=0.0; ioyz[2]=7.5625e-6; ixyz[0]=2.333e-7; ixyz[1]=3.048e-5; ixyz[2]=6.25e-5;
01067
01068
01069
01070
01071 e=d[0][0]; g=d[1][1];
01072
01073
01074 if(beta[0]<Mp->zero) gy=0.0;
01075 else gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
01076
01077 if(beta[1]<Mp->zero) gz=0.;
01078 else gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
01079 gy=0.; gz=0.;
01080 kakr=0.;
01081
01082 eh=pow(g*ixyz[0]/e/ioyz[0],0.5);
01083
01084 eh=eh*dl;
01085 s=sinh(eh); c=cosh(eh);
01086 eh=eh/dl;
01087 ck_matrix (ck, s,c,eh,dl);
01088 xs=0.0;
01089 for (i=0;i<2;i++){
01090
01091 ehx=eh*xs*dl;
01092 s=sinh(ehx); c=cosh(ehx);
01093
01094
01095
01096 bb[0][0]=-e*a/dl;
01097 bb[0][7]= e*a/dl;
01098
01099 ajy= (1.+2.*gy); ajz= (1.+2.*gz); ajk= (1.+2.*kakr);
01100 dderno[0]=(-6.+12.*xs)/ll; dderno[1]=(-2.*(2.+kakr)+6.*xs)/dl; dderno[2]=( 6.-12.*xs)/ll; dderno[3]=(-2.*(1.-kakr)+6.*xs)/dl;
01101
01102 dderny[0]=(-6.+12.*xs)/ll; dderny[1]=(-2.*(2.+gy)+6.*xs)/dl; dderny[2]=( 6.-12.*xs)/ll; dderny[3]=(-2.*(1.-gy)+6.*xs)/dl;
01103
01104
01105 bb[1][2]= -e*ixyz[1]/ajy*dderny[0];
01106 bb[1][3]= -e*ioyz[1]/ajk*dderno[0];
01107 bb[1][4]= e*ixyz[1]/ajy*dderny[1];
01108 bb[1][6]= -e*ioyz[1]/ajk*dderno[1];
01109 bb[1][9]= -e*ixyz[1]/ajy*dderny[2];
01110 bb[1][10]= -e*ioyz[1]/ajk*dderno[2];
01111 bb[1][11]= e*ixyz[1]/ajy*dderny[3];
01112 bb[1][13]= -e*ioyz[1]/ajk*dderno[3];
01113
01114 bb[5][2]= e*ixyz[1]/ajy*12./ll/dl;
01115 bb[5][3]= e*ioyz[1]/ajk*12./ll/dl;
01116 bb[5][4]= -e*ixyz[1]/ajy*6./ll;
01117 bb[5][6]= e*ioyz[1]/ajk*6./ll;
01118 bb[5][9]= e*ixyz[1]/ajy*(-12.)/ll/dl;
01119 bb[5][10]= e*ioyz[1]/ajk*(-12.)/ll/dl;
01120 bb[5][11]=-e*ixyz[1]/ajy*6./ll;
01121 bb[5][13]= e*ioyz[1]/ajk*6./ll;
01122
01123
01124 ddernz[0]=(-6.+12.*xs)/ll; ddernz[1]=(-2.*(2.+gz)+6.*xs)/dl; ddernz[2]=( 6.-12.*xs)/ll; ddernz[3]=(-2.*(1.-gz)+6.*xs)/dl;
01125 bb[2][1]= e*ixyz[2]/ajz*ddernz[0];
01126 bb[2][3]= e*ioyz[2]/ajk*dderno[0];
01127 bb[2][5]= e*ixyz[2]/ajz*ddernz[1];
01128 bb[2][6]= e*ioyz[2]/ajk*dderno[1];
01129 bb[2][ 8]= e*ixyz[2]/ajz*ddernz[2];
01130 bb[2][10]= e*ioyz[2]/ajk*dderno[2];
01131 bb[2][12]= e*ixyz[2]/ajz*ddernz[3];
01132 bb[2][13]= e*ioyz[2]/ajk*dderno[3];
01133
01134 bb[4][1]= e*ixyz[2]/ajz*12./ll/dl;
01135 bb[4][3]= e*ioyz[2]/ajk*12./ll/dl;
01136 bb[4][5]= e*ixyz[2]/ajz*6./ll;
01137 bb[4][6]= e*ioyz[2]/ajk*6./ll;
01138 bb[4][8] = e*ixyz[2]/ajz*(-12.)/ll/dl;
01139 bb[4][10]= e*ioyz[2]/ajk*(-12.)/ll/dl;
01140 bb[4][12]= e*ixyz[2]/ajz*6./ll;
01141 bb[4][13]= e*ioyz[2]/ajk*6./ll;
01142
01143
01144 bb[3][1]= -e*ioyz[2]/ajz*ddernz[0];
01145
01146 bb[3][2]= -e*ioyz[1]/ajy*dderny[0];
01147
01148 bb[3][3]= -e*ioyz[0]*eh*eh*(ck[0][2]*c+ck[0][3]*s);
01149 bb[3][4]= e*ioyz[1]/ajy*dderny[1];
01150
01151 bb[3][5]= -e*ioyz[2]/ajz*ddernz[1];
01152
01153 bb[3][6]= -e*ioyz[0]*eh*eh*(ck[1][2]*c+ck[1][3]*s);
01154 bb[3][8]= -e*ioyz[2]/ajz*ddernz[2];
01155
01156 bb[3][9]= -e*ioyz[1]/ajy*dderny[2];
01157
01158 bb[3][10]=-e*ioyz[0]*eh*eh*(ck[2][2]*c+ck[2][3]*s);
01159 bb[3][11]= e*ioyz[1]/ajy*dderny[3];
01160
01161 bb[3][12]=-e*ioyz[2]/ajz*ddernz[3];
01162
01163 bb[3][13]=-e*ioyz[0]*eh*eh*(ck[3][2]*c+ck[3][3]*s);
01164
01165
01166
01167 bb[6][1]= -e*ioyz[2]/ajz*12./ll/dl;
01168 bb[6][2]= -e*ioyz[1]/ajy*12./ll/dl;
01169 bb[6][3]= -e*ioyz[0]*eh*eh*eh*(ck[0][2]*s+ck[0][3]*c);
01170 bb[6][4]= e*ioyz[1]/ajy*6./ll;
01171 bb[6][5]= -e*ioyz[2]/ajz*6./ll;
01172 bb[6][6]= -e*ioyz[0]*eh*eh*eh*(ck[1][2]*s+ck[1][3]*c);
01173 bb[6][8]= -e*ioyz[2]/ajz*(-12.)/ll/dl;
01174 bb[6][9]= -e*ioyz[1]/ajy*(-12.)/ll/dl;
01175 bb[6][10]=-e*ioyz[0]*eh*eh*eh*(ck[2][2]*s+ck[2][3]*c);
01176 bb[6][11]= e*ioyz[1]/ajy*6./ll;
01177 bb[6][12]=-e*ioyz[2]/ajz*6./ll;
01178 bb[6][13]=-e*ioyz[0]*eh*eh*eh*(ck[3][2]*s+ck[3][3]*c);
01179
01180 bb[7][3]= g*ixyz[0]*(ck[0][1]+eh*(ck[0][2]*s+ck[0][3]*c));
01181 bb[7][6]= g*ixyz[0]*(ck[1][1]+eh*(ck[1][2]*s+ck[1][3]*c));
01182 bb[7][10]= g*ixyz[0]*(ck[2][1]+eh*(ck[2][2]*s+ck[2][3]*c));
01183 bb[7][13]= g*ixyz[0]*(ck[3][1]+eh*(ck[3][2]*s+ck[3][3]*c));
01184
01185
01186 mxv (bb,rl,fx);
01187
01188 xs=1.0;
01189 fprintf (Out,"\n\n MT prvek beam cislo %ld",lcid);
01190 fprintf (Out," \n u, v, w, fx, fy, fz, o \n ");
01191 for (j=0;j<ndofe;j++){
01192 fprintf (Out," %15.5e",rl[j]);
01193 }
01194 fprintf (Out," \n N, My, Mz, B, Vy, Vz, Mo Mt\n ");
01195 for (j=0;j<tncomp;j++){
01196 fprintf (Out," %15.5e",fx[j]);
01197 }
01198 }
01199
01200 }
01201
01202
01203
01204 void beamgen3d::internal_forces1 (long lcid,long eid,long ri,long ci,vector &ifor)
01205 {
01206 long i,j,integr=2,ipp,transf;
01207 double dl,ll,a,ixyz[3],ioyz[3],beta[2],gy,gz,e,g;
01208 double ajy,ajz,ajk, xs, kakr ;
01209 ivector nodes(nne);
01210 vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),f(ndofe),w(integr),gp(integr), ddern(4),dderno(4);
01211 vector fx(4);
01212 matrix n(tncomp,ndofe),d(tncomp,tncomp),tran(3,3),r(tncomp,2);
01213 matrix bb(4,ndofe);
01214
01215
01216
01217 Mt->give_elemnodes (eid,nodes);
01218 Mt->give_node_coord3d (x,y,z,eid);
01219 eldispl (lcid,eid,rl.a);
01220 Mc->give_vectorlcs (eid, vec);
01221 beam_transf_matrix (tran,dl, vec,x,y,z,eid);
01222 ll=dl*dl;
01223
01224
01225 fillm (0.0,r);
01226
01227
01228 transf = Mt->locsystems (nodes);
01229 if (transf>0){
01230 matrix tmat (ndofe,ndofe);
01231 transf_matrix (nodes,tmat);
01232
01233 lgvectortransf (rg,rl,tmat);
01234 }
01235 else{
01236 copyv (rl,rg);
01237 }
01238
01239
01240
01241
01242
01243
01244 ipp=Mt->elements[eid].ipp[0][0];
01245 Mm->matstiff (d,ipp);
01246 Mc->give_area (eid,a);
01247 Mc->give_mominer (eid,ixyz);
01248 Mc->give_shearcoeff (eid,beta);
01249 ioyz[0]=4.2692e-6 ; ioyz[1]=-2.6687e-5; ioyz[2]=0.0; ixyz[0]=2.646e-7; ixyz[1]=2.135e-4; ixyz[2]=3.333e-5;
01250
01251
01252
01253
01254
01255 e=d[0][0]; g=d[1][1];
01256
01257
01258 if(beta[0]<Mp->zero) gy=0.0;
01259 else gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
01260
01261 if(beta[1]<Mp->zero) gz=0.;
01262 else gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
01263
01264 gy=0.; gz=0.; kakr=0.;
01265
01266 xs=0.0;
01267 for (i=0;i<2;i++){
01268
01269 bb[0][0]=-e*a/dl;
01270 bb[0][7]= e*a/dl;
01271
01272 ajy= (1.+2.*gy); ajz= (1.+2.*gz); ajk= (1.+2.*kakr);
01273 dderno[0]=(-6.+12.*xs)/ll; dderno[1]=(-2.*(2.+kakr)+6.*xs)/dl; dderno[2]=( 6.-12.*xs)/ll; dderno[3]=(-2.*(1.-kakr)+6.*xs)/dl;
01274 ddern[0]=(-6.+12.*xs)/ll; ddern[1]=(-2.*(2.+gy)+6.*xs)/dl; ddern[2]=( 6.-12.*xs)/ll; ddern[3]=(-2.*(1.-gy)+6.*xs)/dl;
01275
01276 bb[1][2]= -e*ixyz[1]/ajy*ddern[0];
01277 bb[1][3]= -e*ioyz[1]/ajk*dderno[0];
01278 bb[1][4]= e*ixyz[1]/ajy*ddern[1];
01279 bb[1][6]= -e*ioyz[1]/ajk*dderno[1];
01280 bb[1][9]= -e*ixyz[1]/ajy*ddern[2];
01281 bb[1][10]= -e*ioyz[1]/ajk*dderno[2];
01282 bb[1][11]= e*ixyz[1]/ajy*ddern[3];
01283 bb[1][13]= -e*ioyz[1]/ajk*dderno[3];
01284
01285 bb[3][2]= -e*ioyz[1]/ajy*ddern[0];
01286 bb[3][3]= -e*ioyz[0]/ajk*dderno[0];
01287 bb[3][4]= e*ioyz[1]/ajy*ddern[1];
01288 bb[3][6]= -e*ioyz[0]/ajk*dderno[1];
01289 bb[3][9]= -e*ioyz[1]/ajy*ddern[2];
01290 bb[3][10]= -e*ioyz[0]/ajk*dderno[2];
01291 bb[3][11]= e*ioyz[1]/ajy*ddern[3];
01292 bb[3][13]= -e*ioyz[0]/ajk*dderno[3];
01293
01294
01295 ddern[0]=(-6.+12.*xs)/ll; ddern[1]=(-2.*(2.+gz)+6.*xs)/dl; ddern[2]=( 6.-12.*xs)/ll; ddern[3]=(-2.*(1.-gz)+6.*xs)/dl;
01296 bb[2][1]= e*ixyz[2]/ajz*ddern[0];
01297 bb[2][3]= e*ioyz[2]/ajk*dderno[0];
01298 bb[2][5]= e*ixyz[2]/ajz*ddern[1];
01299 bb[2][6]= e*ioyz[2]/ajk*dderno[1];
01300 bb[2][ 8]= e*ixyz[2]/ajz*ddern[2];
01301 bb[2][10]= e*ioyz[2]/ajk*dderno[2];
01302 bb[2][12]= e*ixyz[2]/ajz*ddern[3];
01303 bb[2][13]= e*ioyz[2]/ajk*dderno[3];
01304
01305 bb[3][1]= -e*ioyz[2]/ajz*ddern[0];
01306 bb[3][5]= -e*ioyz[2]/ajz*ddern[1];
01307 bb[3][8]= -e*ioyz[2]/ajz*ddern[2];
01308 bb[3][12]= -e*ioyz[2]/ajz*ddern[3];
01309
01310 mxv (bb,rl,fx);
01311
01312 xs=1.0;
01313 fprintf (Out,"\n\n MT prvek beam cislo %ld",lcid);
01314 fprintf (Out," \n u, v, w, fx, fy, fz, o \n ");
01315 for (j=0;j<ndofe;j++){
01316 fprintf (Out," %15.5e",rl[j]);
01317 }
01318 fprintf (Out," \n N, My, Mz, B\n ");
01319 for (j=0;j<4;j++){
01320 fprintf (Out," %15.5e",fx[j]);
01321 }
01322 }
01323 }