00001 #include "beam2dspec.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 "normmat.h"
00009 #include <math.h>
00010
00011
00012 beam2dspec::beam2dspec (void)
00013 {
00014 long i;
00015
00016 nne=2; ndofe=6; tncomp=6; napfun=3;
00017 nb=1; intordmm=0; intordism=0; ssst=plbeam;
00018
00019 ncomp = new long [nb];
00020 ncomp[0]=6;
00021
00022 cncomp = new long [nb];
00023 cncomp[0]=0;
00024
00025 nip = new long* [nb];
00026 for (i=0;i<nb;i++){
00027 nip[i] = new long [nb];
00028 }
00029 nip[0][0]=2;
00030 tnip=2;
00031
00032 }
00033
00034 beam2dspec::~beam2dspec (void)
00035 {
00036 long i;
00037
00038 for (i=0;i<nb;i++){
00039 delete [] nip[i];
00040 }
00041 delete [] nip;
00042 }
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 void beam2dspec::transf_matrix (ivector &nodes,matrix &tmat)
00055 {
00056 long i,n,m;
00057
00058 fillm (0.0,tmat);
00059
00060 n=nodes.n;
00061 m=tmat.m;
00062 for (i=0;i<m;i++){
00063 tmat[i][i]=1.0;
00064 }
00065
00066 for (i=0;i<n;i++){
00067 if (Mt->nodes[nodes[i]].transf>0){
00068 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;
00069 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;
00070 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;
00071
00072 }
00073 }
00074 }
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085 void beam2dspec::beam_transf_matrix (double c,double s,matrix &tmat)
00086 {
00087 fillm (0.0,tmat);
00088
00089 tmat[0][0] = c; tmat[0][1] = -1.0*s; tmat[0][2] = 0.0;
00090 tmat[1][0] = s; tmat[1][1] = c; tmat[1][2] = 0.0;
00091 tmat[2][0] = 0.0; tmat[2][1] = 0.0; tmat[2][2] = 1.0;
00092
00093 tmat[3][3] = c; tmat[3][4] = -1.0*s; tmat[3][5] = 0.0;
00094 tmat[4][3] = s; tmat[4][4] = c; tmat[4][5] = 0.0;
00095 tmat[5][3] = 0.0; tmat[5][4] = 0.0; tmat[5][5] = 1.0;
00096 }
00097
00098
00099
00100
00101 void beam2dspec::res_stiffness_matrix (long eid,matrix &sm)
00102 {
00103
00104 stiffness_matrix_expl (eid,0,0,sm);
00105 }
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117 void beam2dspec::stiffness_matrix_expl (long eid,long ri,long ci,matrix &sm)
00118 {
00119 long i,ipp,transf;
00120 double l,s,c;
00121 ivector nodes(nne);
00122 vector x(nne),z(nne);
00123 matrix d(tncomp,tncomp);
00124
00125 Mt->give_elemnodes (eid,nodes);
00126 Mt->give_node_coord2dxz (x,z,eid);
00127
00128
00129 l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
00130 if (l<Mp->zero){
00131 fprintf (stderr,"\n\n zero length of the %ld beam2dspec element",eid);
00132 fprintf (stderr,"\n in function beam2dspec::stiffness_matrix (file %s, line %d).\n",__FILE__,__LINE__);
00133 }
00134
00135 s=(z[1]-z[0])/l;
00136 c=(x[1]-x[0])/l;
00137
00138 fillm (0.0,sm);
00139 ipp=Mt->elements[eid].ipp[ri][ci];
00140
00141
00142 i = Mm->ip[ipp].idm[0];
00143
00144 Mm->normm[i].stiffness_matrix (sm,ipp,l);
00145
00146
00147
00148
00149 matrix tmat (ndofe,ndofe);
00150 beam_transf_matrix (c,s,tmat);
00151 lgmatrixtransf (sm,tmat);
00152
00153 transf = Mt->locsystems (nodes);
00154 if (transf>0){
00155 transf_matrix (nodes,tmat);
00156 glmatrixtransf (sm,tmat);
00157 }
00158
00159 }
00160
00161 void beam2dspec::res_mass_matrix (long eid,matrix &mm)
00162 {
00163 mass_matrix_expl (eid,0,0,mm);
00164 }
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 void beam2dspec::mass_matrix_expl (long eid,long ri,long ci,matrix &mm)
00176 {
00177 long ii,transf;
00178 double e,g,shearcoeff,j,iy,a,l,c,s,kappa,rho;
00179 ivector nodes(nne);
00180 vector x(nne),z(nne);
00181 matrix d (tncomp,tncomp);
00182
00183 Mt->give_elemnodes (eid,nodes);
00184 Mt->give_node_coord2dxz (x,z,eid);
00185
00186 l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
00187 if (l<Mp->zero){
00188 fprintf (stderr,"\n\n zero length of the %ld beam2dspec element",eid);
00189 fprintf (stderr,"\n in function beam2dspec::mass_matrix (file %s, line %d).\n",__FILE__,__LINE__);
00190 }
00191
00192 s=(z[1]-z[0])/l;
00193 c=(x[1]-x[0])/l;
00194
00195 fillm (0.0,mm);
00196 ii=Mt->elements[eid].ipp[ri][ci];
00197 Mc->give_area (eid,a);
00198 Mc->give_shearcoeff (eid,&shearcoeff);
00199 Mc->give_densitye (eid,rho);
00200 Mc->give_mominer (eid,&iy);
00201 Mm->matstiff (d,ii);
00202
00203 e=d[0][0]; g=d[1][1];
00204 if (shearcoeff<Mp->zero) kappa=0.0;
00205 else kappa=6.0*e*iy/shearcoeff/g/a/l/l;
00206
00207 j = (1.0+2.0*kappa)*(1.0+2.0*kappa);
00208
00209 mm[0][0]=rho*a*l/3.0;
00210 mm[0][3]=rho*a*l/6.0;
00211
00212 mm[1][1]=rho*a*l*(13.0/35.0+7.0/5.0*kappa+4.0/3.0*kappa*kappa)/j;
00213 mm[1][2]=rho*a*l*l*(-11.0/210.0-11.0/60*kappa-1.0/6.0*kappa*kappa)/j;
00214 mm[1][4]=rho*a*l*(9.0/70+3.0/5.0*kappa+2.0/3.0*kappa*kappa)/j;
00215 mm[1][5]=rho*a*l*l*(13.0/420+3.0/20.0*kappa+1.0/6.0*kappa*kappa)/j;
00216
00217 mm[2][1]=mm[1][2];
00218 mm[2][2]=rho*a*l*l*l*(1.0/105.0+1.0/30.0*kappa+1.0/30.0*kappa*kappa)/j;
00219 mm[2][4]=rho*a*l*l*(-13.0/420.0-3.0/20.0*kappa-1.0/6.0*kappa*kappa)/j;
00220 mm[2][5]=rho*a*l*l*l*(-1.0/140.0-1.0/30.0*kappa-1.0/30.0*kappa*kappa)/j;
00221
00222 mm[3][0]=mm[0][3];
00223 mm[3][3]=mm[0][0];
00224
00225 mm[4][1]=mm[1][4];
00226 mm[4][2]=mm[2][4];
00227 mm[4][4]=rho*a*l*(13.0/35.0+7.0/5.0*kappa+4.0/3.0*kappa*kappa)/j;
00228 mm[4][5]=rho*a*l*l*(11.0/210.0+11.0/60.0*kappa+1.0/6.0*kappa*kappa)/j;
00229
00230 mm[5][1]=mm[1][5];
00231 mm[5][2]=mm[2][5];
00232 mm[5][4]=mm[4][5];
00233 mm[5][5]=rho*a*l*l*l*(1.0/105+1.0/30.0*kappa+1.0/30.0*kappa*kappa)/j;
00234
00235
00236 matrix tmat (ndofe,ndofe);
00237 beam_transf_matrix (c,s,tmat);
00238 lgmatrixtransf (mm,tmat);
00239
00240 transf = Mt->locsystems (nodes);
00241 if (transf>0){
00242 transf_matrix (nodes,tmat);
00243 glmatrixtransf (mm,tmat);
00244 }
00245 }
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 void beam2dspec::internal_forces (long lcid,long eid,vector &ifor)
00259 {
00260 long i,ipp;
00261 double l,c,s;
00262 vector x(nne),z(nne),lifor(ndofe),r(ndofe);
00263 matrix d(tncomp,tncomp);
00264
00265
00266 Mt->give_node_coord2dxz (x,z,eid);
00267
00268
00269 eldispl (lcid,eid,r.a);
00270
00271
00272 l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
00273 if (l<Mp->zero){
00274 fprintf (stderr,"\n\n zero length of the %ld beam2dspec element",eid);
00275 fprintf (stderr,"\n in function beam2dspec::internal_forces (file %s, line %d).\n",__FILE__,__LINE__);
00276 }
00277
00278 s=(z[1]-z[0])/l;
00279 c=(x[1]-x[0])/l;
00280
00281
00282 fillv (0.0,ifor);
00283
00284
00285 ipp=Mt->elements[eid].ipp[0][0];
00286
00287 i=Mm->ip[ipp].idm[0];
00288
00289
00290
00291 Mm->normm[i].internal_forces (ifor,ipp,l);
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301 }
00302
00303 void beam2dspec::res_internal_forces (long lcid,long eid,vector &ifor)
00304 {
00305 internal_forces (lcid,eid,ifor);
00306 }
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318 void beam2dspec::stresses (long eid,long lcid)
00319 {
00320 long i,ipp;
00321 vector x(nne),z(nne),f(ndofe),r(ndofe);
00322 matrix d(tncomp,tncomp);
00323
00324
00325 Mt->give_node_coord2dxz (x,z,eid);
00326
00327
00328 eldispl (lcid,eid,r.a);
00329
00330
00331 ipp=Mt->elements[eid].ipp[0][0];
00332
00333 i=Mm->ip[ipp].idm[0];
00334
00335 matrix sm(6,6);
00336 stiffness_matrix_expl (eid,0,0,sm);
00337
00338 mxv (sm,r,f);
00339
00340 for (i=0;i<6;i++){
00341 Mm->ip[ipp].stress[i]=f[i];
00342 }
00343
00344 }
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 void beam2dspec::strains (long eid,long lcid)
00355 {
00356 long ii,transf;
00357 double l,s,c;
00358 ivector nodes(nne);
00359 vector x(nne),z(nne),rl(ndofe),rg(ndofe);
00360
00361 Mt->give_elemnodes (eid,nodes);
00362 Mt->give_node_coord2dxz (x,z,eid);
00363
00364
00365 l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
00366 if (l<Mp->zero){
00367 fprintf (stderr,"\n\n zero length of the %ld beamel2d element",eid);
00368 fprintf (stderr,"\n in function beamel2d::strains (file %s, line %d).\n",__FILE__,__LINE__);
00369 }
00370
00371 s=(z[1]-z[0])/l;
00372 c=(x[1]-x[0])/l;
00373
00374 eldispl (lcid,eid,rl.a);
00375
00376
00377 matrix tmat (ndofe,ndofe);
00378 transf = Mt->locsystems (nodes);
00379 if (transf>0){
00380 transf_matrix (nodes,tmat);
00381
00382 lgvectortransf (rg,rl,tmat);
00383 }
00384 else{
00385 copyv (rl,rg);
00386 }
00387 beam_transf_matrix (c,s,tmat);
00388
00389 glvectortransf (rg,rl,tmat);
00390
00391
00392 ii=Mt->elements[eid].ipp[0][0];
00393 Mm->storestrain (lcid,ii,rl);
00394 }
00395
00396
00397
00398
00399
00400
00401
00402 void beam2dspec::res_mainip_strains (long lcid,long eid)
00403 {
00404 long ipp;
00405 vector aux,r(ndofe),eps(3);
00406 ivector nodes(nne);
00407 matrix tmat;
00408
00409
00410 eldispl (lcid,eid,r.a);
00411
00412
00413 long transf = Mt->locsystems (nodes);
00414 if (transf>0){
00415 allocv (ndofe,aux);
00416 allocm (ndofe,ndofe,tmat);
00417 transf_matrix (nodes,tmat);
00418
00419 lgvectortransf (aux,r,tmat);
00420 copyv (aux,r);
00421 destrv (aux);
00422 destrm (tmat);
00423 }
00424
00425 eps[0]=r[0];
00426 eps[1]=r[1];
00427 eps[2]=r[2];
00428
00429
00430 ipp=Mt->elements[eid].ipp[0][0];
00431
00432 Mm->storestrain (lcid,ipp,0,eps);
00433
00434 }
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444 void beam2dspec::res_mainip_stresses (long lcid,long eid)
00445 {
00446 double s,c,l;
00447 ivector nodes(nne);
00448 vector x(nne),z(nne),r(ndofe),gf(ndofe),lf(ndofe),aux;
00449 matrix sm(ndofe,ndofe),tmat(ndofe,ndofe);
00450
00451 eldispl (lcid,eid,r.a);
00452
00453
00454 Mt->give_elemnodes (eid,nodes);
00455 long transf = Mt->locsystems (nodes);
00456 if (transf>0){
00457 allocv (ndofe,aux);
00458
00459 transf_matrix (nodes,tmat);
00460
00461 lgvectortransf (aux,r,tmat);
00462 copyv (aux,r);
00463 destrv (aux);
00464
00465 }
00466
00467
00468
00469 stiffness_matrix_expl (eid,0,0,sm);
00470
00471 mxv (sm,r,gf);
00472
00473
00474
00475
00476 Mt->give_node_coord2dxz (x,z,eid);
00477
00478
00479 l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
00480 if (l<Mp->zero){
00481 fprintf (stderr,"\n\n zero length of the %ld beamel2d element",eid);
00482 fprintf (stderr,"\n in function beamel2d::res_mainip_stresses (file %s, line %d).\n",__FILE__,__LINE__);
00483 }
00484
00485
00486 s=(z[1]-z[0])/l;
00487 c=(x[1]-x[0])/l;
00488
00489
00490
00491
00492 beam_transf_matrix (c,s,tmat);
00493
00494
00495
00496 lgvectortransf (gf,lf,tmat);
00497
00498 long ipp;
00499 ipp=Mt->elements[eid].ipp[0][0];
00500 Mm->ip[ipp].stress[0]=lf[0]*(-1.0);
00501 Mm->ip[ipp].stress[1]=lf[1]*(-1.0);
00502 Mm->ip[ipp].stress[2]=lf[2]*(-1.0);
00503
00504 ipp++;
00505
00506
00507 Mm->ip[ipp].stress[0]=lf[3];
00508 Mm->ip[ipp].stress[1]=lf[4];
00509 Mm->ip[ipp].stress[2]=lf[5];
00510
00511
00512 }