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