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