00001 #include <math.h>
00002 #include "barel3d.h"
00003 #include "global.h"
00004 #include "globmat.h"
00005 #include "genfile.h"
00006 #include "element.h"
00007 #include "node.h"
00008 #include "intpoints.h"
00009 #include "loadcase.h"
00010
00011
00012
00013
00014
00015
00016 barel3d::barel3d (void)
00017 {
00018 long i;
00019
00020
00021 nne=2;
00022
00023 ndofe=6;
00024
00025 tncomp=1;
00026
00027 napfun=3;
00028 ssst = bar;
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 nip[0][0]=1;
00051
00052
00053 tnip=0;
00054 for (i=0;i<nb;i++){
00055 tnip+=nip[i][i];
00056 }
00057
00058 intordsm[0][0]=1;
00059 }
00060
00061
00062
00063 barel3d::~barel3d (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 void barel3d::dirvect (vector &s,double &l,vector &x,vector &y,vector &z)
00093 {
00094 s[0]=x[1]-x[0];
00095 s[1]=y[1]-y[0];
00096 s[2]=z[1]-z[0];
00097
00098 l=sqrt(s[0]*s[0]+s[1]*s[1]+s[2]*s[2]);
00099
00100 if (l<Mp->zero){
00101 fprintf (stderr,"\n\n zero norm of direction vector in function dirvect (file %s, line %d).\n",__FILE__,__LINE__);
00102 }
00103
00104 s[0]/=l;
00105 s[1]/=l;
00106 s[2]/=l;
00107 }
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 void barel3d::geom_matrix (matrix &gm,vector &s,double l)
00123 {
00124 gm[0][0]=s[0]/l*(-1.0);
00125 gm[0][1]=s[1]/l*(-1.0);
00126 gm[0][2]=s[2]/l*(-1.0);
00127 gm[0][3]=s[0]/l;
00128 gm[0][4]=s[1]/l;
00129 gm[0][5]=s[2]/l;
00130 }
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 double barel3d::approx (double xi,vector &nodval)
00145 {
00146 double f;
00147 vector bf(nne);
00148
00149 bf_lin_1d (bf.a,xi);
00150
00151 scprd (bf,nodval,f);
00152
00153 return f;
00154 }
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 void barel3d::transf_matrix (ivector &nodes,matrix &tmat)
00170 {
00171 long i;
00172 fillm (0.0,tmat);
00173
00174 for (i=0;i<tmat.m;i++){
00175 tmat[i][i]=1.0;
00176 }
00177
00178 for (i=0;i<nodes.n;i++){
00179 if (Mt->nodes[nodes[i]].transf>0){
00180 tmat[i*3+0][i*3]=Mt->nodes[nodes[i]].e1[0];
00181 tmat[i*3+1][i*3]=Mt->nodes[nodes[i]].e1[1];
00182 tmat[i*3+2][i*3]=Mt->nodes[nodes[i]].e1[2];
00183
00184 tmat[i*3+0][i*3+1]=Mt->nodes[nodes[i]].e2[0];
00185 tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e2[1];
00186 tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e2[2];
00187
00188 tmat[i*3+0][i*3+2]=Mt->nodes[nodes[i]].e3[0];
00189 tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e3[1];
00190 tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e3[2];
00191 }
00192 }
00193 }
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211 void barel3d::tran_mat (vector &x, matrix &tmat,vector &gx,vector &gy)
00212 {
00213 long i,j;
00214 vector gv(2*gx.n),lv(2*gx.n);
00215
00216 j=0;
00217 for (i=0;i<gx.n;i++){
00218 gv[j]=gx[i]; j++;
00219 gv[j]=gy[i]; j++;
00220 }
00221
00222
00223 glvectortransf (gv,lv,tmat);
00224
00225 for (i=0;i<gx.n;i++){
00226 x[i]=lv[i*2];
00227 }
00228 }
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 void barel3d::stiffness_matrix (long eid,long ri,long ci,vector &s,matrix &sm)
00242 {
00243 long ipp;
00244 double l,a,jac;
00245 vector x(nne),y(nne),z(nne);
00246 matrix gm(tncomp,ndofe),d(tncomp,tncomp);
00247
00248 fillm (0.0,sm);
00249
00250
00251 Mt->give_node_coord3d (x,y,z,eid);
00252
00253
00254 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]));
00255 if (l<Mp->zero)
00256 print_err("zero length of the truss element",__FILE__,__LINE__,__func__);
00257
00258
00259
00260 geom_matrix (gm,s,l);
00261
00262 ipp=Mt->elements[eid].ipp[ri][ci];
00263
00264 Mm->matstiff (d,ipp);
00265
00266 Mc->give_area (eid,a);
00267 jac=a*l;
00268
00269 bdbj (sm.a,gm.a,d.a,jac,gm.m,gm.n);
00270 }
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 void barel3d::res_stiffness_matrix (long eid,matrix &sm)
00286 {
00287 long transf;
00288 double l;
00289 ivector nodes(nne);
00290 vector x(nne),y(nne),z(nne),s(3);
00291 matrix tmat (ndofe,ndofe);
00292
00293
00294 Mt->give_node_coord3d (x,y,z,eid);
00295
00296
00297 dirvect (s,l,x,y,z);
00298
00299
00300 stiffness_matrix (eid,0,0,s,sm);
00301
00302
00303
00304 Mt->give_elemnodes (eid,nodes);
00305
00306
00307 transf = Mt->locsystems (nodes);
00308 if (transf>0){
00309 transf_matrix (nodes,tmat);
00310 glmatrixtransf (sm,tmat);
00311 }
00312 }
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326 void barel3d::mass_matrix (long eid,matrix &mm)
00327 {
00328 double l,a,rho;
00329 vector x(nne),y(nne),z(nne);
00330
00331 fillm (0.0,mm);
00332
00333 Mt->give_node_coord3d (x,y,z,eid);
00334
00335
00336 l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
00337 if (l<Mp->zero){
00338 fprintf (stderr,"\n\n zero length of the %ld truss element",eid);
00339 fprintf (stderr,"\n in function barel3d::mass_matrix (%s, line %d).\n",__FILE__,__LINE__);
00340 }
00341 Mc->give_area (eid,a);
00342 Mc->give_densitye (eid,rho);
00343
00344 mm[0][0]=a*rho*l*1.0/3.0;
00345 mm[0][1]=0.0;
00346 mm[0][2]=a*rho*l*1.0/6.0;
00347 mm[0][3]=0.0;
00348
00349 mm[1][0]=0.0;
00350 mm[1][1]=a*rho*l*1.0/3.0;
00351 mm[1][2]=0.0;
00352 mm[1][3]=a*rho*l*1.0/6.0;
00353
00354 mm[2][0]=a*rho*l*1.0/6.0;
00355 mm[2][1]=0.0;
00356 mm[2][2]=a*rho*l*1.0/3.0;
00357 mm[2][3]=0.0;
00358
00359 mm[3][0]=0.0;
00360 mm[3][1]=a*rho*l*1.0/6.0;
00361 mm[3][2]=0.0;
00362 mm[3][3]=a*rho*l*1.0/3.0;
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 void barel3d::res_mass_matrix (long eid,matrix &mm)
00379 {
00380 vector gx(nne),gy(nne),x(nne);
00381 long transf;
00382
00383 Mt->give_node_coord2d (gx,gy,eid);
00384 giveloccoord (gx,gy,x);
00385
00386 mass_matrix (eid,mm);
00387
00388
00389
00390 ivector nodes (nne);
00391 Mt->give_elemnodes (eid,nodes);
00392 transf = Mt->locsystems (nodes);
00393 if (transf>0)
00394 {
00395 matrix tmat (ndofe,ndofe);
00396 transf_matrix (nodes,tmat);
00397 glmatrixtransf (mm,tmat);
00398 }
00399 }
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 void barel3d::giveloccoord(vector &x, vector &y,vector &lx)
00413 {
00414 double l;
00415
00416 lx[0] = 0.0;
00417 l = sqr(x[1] - x[0]) + sqr(y[1] - y[0]);
00418 l = sqrt(l);
00419 lx[1] = l;
00420 }
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 void barel3d::res_ip_strains (long lcid,long eid)
00433 {
00434 long transf;
00435 double l;
00436 vector gx(nne),gy(nne),gz(nne),s(3),r(ndofe),aux(ndofe);
00437 ivector cn(ndofe),nodes(nne);
00438 matrix tmat(ndofe,ndofe);
00439
00440
00441 Mt->give_node_coord3d (gx,gy,gz,eid);
00442
00443 dirvect (s,l,gx,gy,gz);
00444
00445
00446 eldispl (lcid,eid,r.a);
00447
00448
00449
00450
00451 Mt->give_elemnodes (eid,nodes);
00452 transf = Mt->locsystems (nodes);
00453 if (transf>0){
00454
00455 transf_matrix (nodes,tmat);
00456 lgvectortransf (aux,r,tmat);
00457 copyv (aux,r);
00458 }
00459
00460 ip_strains (lcid,eid,0,0,s,r,l);
00461 }
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478 void barel3d::ip_strains (long lcid,long eid,long ri,long ci,vector &s,vector &r,double l)
00479 {
00480 long i,ipid;
00481 vector eps(tncomp);
00482 matrix gm(tncomp,ndofe);
00483
00484 ipid=Mt->elements[eid].ipp[ri][ci];
00485 for (i=0;i<intordsm[0][0];i++){
00486
00487
00488 geom_matrix (gm,s,l);
00489
00490 mxv (gm,r,eps);
00491
00492
00493 Mm->storestrain (lcid,ipid,eps);
00494 ipid++;
00495 }
00496 }
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509 void barel3d::nod_strains_ip (long lcid,long eid,long ri,long ci)
00510 {
00511 long i,j,ipid;
00512 ivector ipnum(nne),nod(nne);
00513 vector eps(tncomp);
00514
00515
00516
00517 ipid=Mt->elements[eid].ipp[ri][ci];
00518 nodip_bar (ipid,intordsm[0][0],ipnum);
00519
00520
00521 Mt->give_elemnodes (eid,nod);
00522
00523 for (i=0;i<nne;i++){
00524
00525 Mm->givestrain (lcid,ipnum[i],eps);
00526
00527
00528 j=nod[i];
00529 Mt->nodes[j].storestrain (lcid,0,eps);
00530 }
00531 }
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544 void barel3d::strains (long lcid,long eid,long ri,long ci)
00545 {
00546 long i,naep,ncp,sid;
00547 vector coord,eps;
00548
00549 switch (Mm->stra.tape[eid]){
00550 case nowhere:{
00551 break;
00552 }
00553 case intpts:{
00554
00555 break;
00556 }
00557 case enodes:{
00558 nod_strains_ip (lcid,eid,ri,ci);
00559 break;
00560 }
00561 case userdefined:{
00562
00563 naep = Mm->stra.give_naep (eid);
00564 ncp = Mm->stra.give_ncomp (eid);
00565 sid = Mm->stra.give_sid (eid);
00566 allocv (ncp,eps);
00567 allocv (1,coord);
00568 for (i=0;i<naep;i++){
00569 Mm->stra.give_aepcoord (sid,i,coord);
00570
00571 if (Mp->strainaver==0)
00572
00573 if (Mp->strainaver==1)
00574
00575
00576 Mm->stra.storevalues(lcid,eid,i,eps);
00577 }
00578 destrv (eps);
00579 destrv (coord);
00580 break;
00581 }
00582 default:
00583 print_err("unknown strain point is required",__FILE__,__LINE__,__func__);
00584 }
00585
00586 }
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 void barel3d::res_ip_stresses (long lcid,long eid)
00599 {
00600 ip_stresses (lcid,eid,0,0);
00601 }
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 void barel3d::ip_stresses (long lcid,long eid,long ri,long ci)
00615 {
00616 long ipid;
00617
00618 ipid=Mt->elements[eid].ipp[ri][ci];
00619
00620
00621 if (Mp->strcomp==1)
00622 Mm->computenlstresses (ipid);
00623 }
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636 void barel3d::ip_elast_stresses (long lcid,long eid,long ri,long ci)
00637 {
00638 long i,ipid;
00639 vector eps(tncomp),sig(tncomp);
00640 matrix d(tncomp,tncomp);
00641
00642 ipid=Mt->elements[eid].ipp[ri][ci];
00643
00644 for (i=0;i<intordsm[0][0];i++){
00645
00646 Mm->matstiff (d,ipid);
00647
00648 Mm->givestrain (lcid,ipid,eps);
00649
00650 mxv (d,eps,sig);
00651
00652 Mm->storestress (lcid,ipid,sig);
00653 ipid++;
00654 }
00655 }
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668 void barel3d::nod_stresses_ip (long lcid,long eid,long ri,long ci)
00669 {
00670 long i,j,ipid;
00671 ivector ipnum(nne),nod(nne);
00672 vector sig(tncomp);
00673
00674
00675
00676 ipid=Mt->elements[eid].ipp[ri][ci];
00677 nodip_bar (ipid,intordsm[0][0],ipnum);
00678
00679
00680 Mt->give_elemnodes (eid,nod);
00681
00682 for (i=0;i<nne;i++){
00683
00684 Mm->givestress (lcid,ipnum[i],sig);
00685
00686
00687 j=nod[i];
00688 Mt->nodes[j].storestress (lcid,0,sig);
00689 }
00690 }
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702 void barel3d::stresses (long lcid,long eid,long ri,long ci)
00703 {
00704 long i,naep,ncp,sid;
00705 vector coord,sig;
00706
00707 switch (Mm->stre.tape[eid]){
00708 case nowhere:{
00709 break;
00710 }
00711 case intpts:{
00712
00713 break;
00714 }
00715 case enodes:{
00716 nod_stresses_ip (lcid,eid,ri,ci);
00717 break;
00718 }
00719 case userdefined:{
00720
00721 naep = Mm->stre.give_naep (eid);
00722 ncp = Mm->stre.give_ncomp (eid);
00723 sid = Mm->stre.give_sid (eid);
00724 allocv (ncp,sig);
00725 allocv (1,coord);
00726 for (i=0;i<naep;i++){
00727 Mm->stre.give_aepcoord (sid,i,coord);
00728
00729 if (Mp->stressaver==0)
00730
00731 if (Mp->stressaver==1)
00732
00733
00734 Mm->stre.storevalues(lcid,eid,i,sig);
00735 }
00736 destrv (sig);
00737 destrv (coord);
00738 break;
00739 }
00740 default:
00741 print_err("unknown stress point is required",__FILE__,__LINE__,__func__);
00742 }
00743 }
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755 void barel3d::nod_eqother_ip (long lcid,long eid,long ri,long ci)
00756 {
00757 long i,j,ipid,ncompo;
00758 ivector ipnum(nne),nod(nne);
00759 vector eqother;
00760
00761
00762
00763 ipid=Mt->elements[eid].ipp[ri][ci];
00764 nodip_bar (ipid,intordsm[0][0],ipnum);
00765
00766
00767 Mt->give_elemnodes (eid,nod);
00768
00769 for (i=0;i<nne;i++){
00770 ncompo = Mm->givencompeqother (ipnum[i],0);
00771 allocv (ncompo,eqother);
00772 Mm->giveeqother (ipnum[i],0,ncompo,eqother.a);
00773
00774
00775 j=nod[i];
00776 Mt->nodes[j].storeother (lcid,0,ncompo,eqother);
00777
00778 destrv (eqother);
00779 }
00780 }
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797 void barel3d::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
00798 {
00799 long ii,jj,i,aipp;
00800 double xi;
00801 vector w,gp;
00802 vector x(nne),y(nne),z(nne);
00803
00804 Mt->give_node_coord3d (x,y,z,eid);
00805 for (ii = 0; ii < nb; ii++){
00806 for (jj = 0; jj < nb; jj++){
00807
00808 aipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00809 if (intordsm[ii][jj] == 0)
00810 continue;
00811 allocv (intordsm[ii][jj],gp);
00812 allocv (intordsm[ii][jj],w);
00813 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00814
00815 for (i=0;i<intordsm[ii][jj];i++){
00816 xi=gp[i];
00817 if (aipp==ipp){
00818 coord[0]=approx (xi,x);
00819 coord[1]=approx (xi,y);
00820 coord[2]=approx (xi,z);
00821 }
00822 aipp++;
00823 }
00824 destrv (gp); destrv (w);
00825 }
00826 }
00827 }
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840 void barel3d::ipvolume (long eid,long ri,long ci)
00841 {
00842 long ipp;
00843 double l,a,vol;
00844 vector x(nne),y(nne),z(nne);
00845
00846
00847 Mt->give_node_coord3d (x,y,z,eid);
00848
00849
00850 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]));
00851 if (l<Mp->zero){
00852 char emsg[100];
00853 sprintf(emsg, "zero length of the %ld truss element", eid);
00854 print_err(emsg, __FILE__, __LINE__, __func__);
00855 }
00856
00857 Mc->give_area (eid,a);
00858
00859 vol=a*l;
00860
00861 ipp=Mt->elements[eid].ipp[ri][ci];
00862 Mm->storeipvol (ipp,vol);
00863 }
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878 void barel3d::intpointval (long eid,vector &nodval,vector &ipval)
00879 {
00880 long ii,jj,i,k,ipp;
00881 double xi;
00882 vector w,gp;
00883
00884 k=0;
00885 for (ii = 0; ii < nb; ii++){
00886 for (jj = 0; jj < nb; jj++){
00887
00888
00889 ipp=Mt->elements[eid].ipp[ii][jj];
00890 if (intordsm[ii][jj] == 0)
00891 continue;
00892 allocv (intordsm[ii][jj],gp);
00893 allocv (intordsm[ii][jj],w);
00894 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00895
00896 for (i=0;i<intordsm[ii][jj];i++){
00897 xi=gp[i];
00898
00899 ipval[k] = approx (xi,nodval);
00900 k++;
00901 ipp++;
00902 }
00903 destrv (gp); destrv (w);
00904 }
00905 }
00906 }
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932 void barel3d::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
00933 {
00934 long i, j, k, ipp;
00935 long ii, jj, nv = nodval.n;
00936 double xi, ipval;
00937 vector w, gp, anv(nne);
00938 long nstra, nstre, ncompstr, ncompeqother;
00939 long idstra, idstre, idoth, idic;
00940 inictype ict;
00941
00942 nstra = idstra = nstre = idstre = idoth = idic = 0;
00943
00944 ict = ictn[0];
00945 for (i=0; i<nne; i++)
00946 {
00947 if (ictn[i] != ict)
00948 {
00949 print_err("Incompatible types of initial conditions on element %ld\n"
00950 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
00951 abort();
00952 }
00953 }
00954 for (j = 0; j < nv; j++)
00955 {
00956 for(i = 0; i < nne; i++)
00957 anv[i] = nodval[i][j];
00958 for (ii = 0; ii < nb; ii++)
00959 {
00960 for (jj = 0; jj < nb; jj++)
00961 {
00962 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00963 if (intordsm[ii][jj] == 0)
00964 continue;
00965 allocv (intordsm[ii][jj],gp);
00966 allocv (intordsm[ii][jj],w);
00967 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00968 for (k = 0; k < intordsm[ii][jj]; k++)
00969 {
00970 xi=gp[k];
00971
00972 ipval = approx (xi,anv);
00973 if ((ictn[i] & inistrain) && (j < Mm->ip[ipp].ncompstr))
00974 {
00975 Mm->ip[ipp].strain[idstra] += ipval;
00976 ipp++;
00977 continue;
00978 }
00979 if ((ictn[i] & inistress) && (j < nstra + Mm->ip[ipp].ncompstr))
00980 {
00981 Mm->ip[ipp].stress[idstre] += ipval;
00982 ipp++;
00983 continue;
00984 }
00985 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
00986 {
00987 Mm->ip[ipp].eqother[idoth] += ipval;
00988 ipp++;
00989 continue;
00990 }
00991 if ((ictn[0] & inicond) && (j < nv))
00992 {
00993 if (Mm->ic[ipp] == NULL)
00994 {
00995 Mm->ic[ipp] = new double[nv-j];
00996 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
00997 }
00998 Mm->ic[ipp][idic] += ipval;
00999 ipp++;
01000 continue;
01001 }
01002 ipp++;
01003 }
01004 destrv (gp); destrv (w);
01005 }
01006 }
01007 ipp=Mt->elements[eid].ipp[ri][ci];
01008 ncompstr = Mm->ip[ipp].ncompstr;
01009 ncompeqother = Mm->ip[ipp].ncompeqother;
01010 if ((ictn[0] & inistrain) && (j < ncompstr))
01011 {
01012 nstra++;
01013 idstra++;
01014 continue;
01015 }
01016 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01017 {
01018 nstre++;
01019 idstre++;
01020 continue;
01021 }
01022 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
01023 {
01024 idoth++;
01025 continue;
01026 }
01027 if ((ictn[0] & inicond) && (j < nv))
01028 {
01029 idic++;
01030 continue;
01031 }
01032 }
01033 }
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050 void barel3d::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
01051 {
01052 integratedquant iq;
01053 iq=locstress;
01054 vector x(nne),y(nne),z(nne);
01055
01056
01057 compute_nlstress (lcid,eid,ri,ci);
01058
01059 Mt->give_node_coord3d (x,y,z,eid);
01060
01061
01062 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01063 }
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080 void barel3d::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
01081 {
01082 integratedquant iq=nonlocstress;
01083 vector x(nne),y(nne),z(nne);
01084 Mt->give_node_coord3d (x,y,z,eid);
01085
01086
01087 compute_nonloc_nlstress (lcid,eid,ri,ci);
01088
01089
01090 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01091 }
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107 void barel3d::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
01108 {
01109 integratedquant iq;
01110 iq=stressincr;
01111 vector x(nne),y(nne),z(nne);
01112
01113
01114 compute_nlstressincr (lcid,eid,ri,ci);
01115
01116 Mt->give_node_coord3d (x,y,z,eid);
01117
01118
01119 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01120
01121 }
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137 void barel3d::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor)
01138 {
01139 integratedquant iq;
01140 iq=eigstress;
01141 vector x(nne),y(nne),z(nne);
01142
01143 Mt->give_node_coord3d (x,y,z,eid);
01144
01145
01146 if (Mp->eigstrcomp)
01147 compute_eigstress (lcid,eid,ri,ci);
01148
01149
01150 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y,z);
01151 }
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167 void barel3d::res_internal_forces (long lcid,long eid,vector &ifor)
01168 {
01169 long transf;
01170 ivector nodes(nne);
01171
01172 internal_forces (lcid,eid,0,0,ifor);
01173
01174
01175 Mt->give_elemnodes (eid,nodes);
01176
01177
01178 transf = Mt->locsystems (nodes);
01179 if (transf>0){
01180 matrix tmat (ndofe,ndofe);
01181 vector v(ndofe);
01182 transf_matrix (nodes,tmat);
01183 glvectortransf (ifor,v,tmat);
01184 copyv (v,ifor);
01185 }
01186 }
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201 void barel3d::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
01202 {
01203 long transf;
01204 ivector nodes(nne);
01205
01206 nonloc_internal_forces (lcid,eid,0,0,ifor);
01207
01208
01209 Mt->give_elemnodes (eid,nodes);
01210
01211
01212 transf = Mt->locsystems (nodes);
01213 if (transf>0){
01214 matrix tmat (ndofe,ndofe);
01215 vector v(ndofe);
01216 transf_matrix (nodes,tmat);
01217 glvectortransf (ifor,v,tmat);
01218 copyv (v,ifor);
01219 }
01220 }
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235 void barel3d::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01236 {
01237 long transf;
01238 ivector nodes(nne);
01239
01240 incr_internal_forces (lcid,eid,0,0,ifor);
01241
01242
01243 Mt->give_elemnodes (eid,nodes);
01244
01245
01246 transf = Mt->locsystems (nodes);
01247 if (transf>0){
01248 matrix tmat (ndofe,ndofe);
01249 vector v(ndofe);
01250 transf_matrix (nodes,tmat);
01251 glvectortransf (ifor,v,tmat);
01252 copyv (v,ifor);
01253 }
01254 }
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269 void barel3d::res_eigstrain_forces (long lcid,long eid,vector &nfor)
01270 {
01271 long transf;
01272 ivector nodes(nne);
01273
01274 eigstrain_forces (lcid,eid,0,0,nfor);
01275
01276
01277 Mt->give_elemnodes (eid,nodes);
01278
01279
01280 transf = Mt->locsystems (nodes);
01281 if (transf>0){
01282 matrix tmat (ndofe,ndofe);
01283 vector v(ndofe);
01284 transf_matrix (nodes,tmat);
01285 glvectortransf (nfor,v,tmat);
01286 copyv (v,nfor);
01287 }
01288 }
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301 void barel3d::compute_nlstress (long lcid,long eid,long ri,long ci)
01302 {
01303 long ipid;
01304
01305 ipid=Mt->elements[eid].ipp[ri][ci];
01306
01307
01308 if (Mp->strcomp==1)
01309 Mm->computenlstresses (ipid);
01310
01311 }
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325 void barel3d::local_values (long lcid,long eid,long ri,long ci)
01326 {
01327 long ipp;
01328
01329 ipp=Mt->elements[eid].ipp[ri][ci];
01330
01331
01332
01333 if (Mp->strcomp==1)
01334 Mm->computenlstresses (ipp);
01335
01336 }
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349 void barel3d::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01350 {
01351 long ipid;
01352
01353 ipid=Mt->elements[eid].ipp[ri][ci];
01354
01355
01356 if (Mp->strcomp==1)
01357 Mm->compnonloc_nlstresses (ipid);
01358
01359 }
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372 void barel3d::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01373 {
01374 long ipid;
01375
01376 ipid=Mt->elements[eid].ipp[ri][ci];
01377
01378
01379 if (Mp->strcomp==1)
01380 Mm->computenlstressesincr (ipid);
01381
01382 }
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395 void barel3d::compute_eigstress (long lcid,long eid,long ri,long ci)
01396 {
01397 long ipid,rmid;
01398 vector eigstr(tncomp),sig(tncomp);
01399 matrix d(tncomp,tncomp);
01400
01401
01402 ipid=Mt->elements[eid].ipp[ri][ci];
01403
01404 rmid = Mm->ip[ipid].grmid ();
01405
01406 if (rmid>-1){
01407
01408 Mm->eigenstresses (sig,ipid,rmid);
01409 }else{
01410
01411 Mm->giveeigstrain (ipid,eigstr);
01412
01413
01414 Mm->matstiff (d,ipid);
01415
01416 mxv (d,eigstr,sig);
01417 }
01418
01419 Mm->storeeigstress (ipid,sig);
01420
01421 }
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437 void barel3d::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y,vector &z)
01438 {
01439 long ipp;
01440 double l,a;
01441 ivector nodes(nne);
01442 vector ipv(tncomp),contr(ndofe);
01443 vector s(3);
01444 matrix gm(tncomp,ndofe);
01445
01446
01447 dirvect (s,l,x,y,z);
01448
01449 ipp=Mt->elements[eid].ipp[ri][ci];
01450
01451
01452 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]));
01453 if (l<Mp->zero)
01454 print_err("zero length of the truss element",__FILE__,__LINE__,__func__);
01455
01456
01457 Mm->givequantity (iq,lcid,ipp,cncomp[0],ipv);
01458
01459 geom_matrix (gm,s,l);
01460
01461 mtxv (gm,ipv,contr);
01462
01463 Mc->give_area (eid,a);
01464 cmulv (a*l,contr);
01465
01466 addv(contr,nv,nv);
01467 }
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481 void barel3d::find_extreme_strains (vector &min,vector &max,long lcid,long eid)
01482 {
01483 long i,ipp;
01484 vector eps(tncomp);
01485
01486
01487 ipp=Mt->elements[eid].ipp[0][0];
01488
01489
01490 Mm->givestrain (lcid,ipp,eps);
01491
01492
01493 for (i=0;i<tncomp;i++){
01494 if (min[i]>eps[i])
01495 min[i]=eps[i];
01496 if (max[i]<eps[i])
01497 max[i]=eps[i];
01498 }
01499
01500 }
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512 void barel3d::find_extreme_stresses (vector &min,vector &max,long lcid,long eid)
01513 {
01514 long i,ipp;
01515 vector sig(tncomp);
01516
01517
01518 ipp=Mt->elements[eid].ipp[0][0];
01519
01520
01521 Mm->givestress (lcid,ipp,sig);
01522
01523
01524 for (i=0;i<tncomp;i++){
01525 if (min[i]>sig[i])
01526 min[i]=sig[i];
01527 if (max[i]<sig[i])
01528 max[i]=sig[i];
01529 }
01530
01531 }