00001 #include <math.h>
00002 #include "barelq3d.h"
00003 #include "barel2d.h"
00004 #include "global.h"
00005 #include "globmat.h"
00006 #include "genfile.h"
00007 #include "element.h"
00008 #include "node.h"
00009 #include "intpoints.h"
00010 #include "loadcase.h"
00011
00012
00013
00014 barelq3d::barelq3d (void)
00015 {
00016 long i;
00017
00018 nne=3; ndofe=9; tncomp=1; napfun=1;
00019 tnip=3;
00020 ssst = bar;
00021
00022 nb=1;
00023
00024 ncomp = new long [nb];
00025 ncomp[0]=1;
00026
00027 cncomp = new long [nb];
00028 cncomp[0]=0;
00029
00030 nip = new long* [nb];
00031 intordsm = new long* [nb];
00032 for (i=0;i<nb;i++){
00033 nip[i] = new long [nb];
00034 intordsm[i] = new long [nb];
00035 }
00036
00037 nip[0][0]=3;
00038
00039 intordsm[0][0]=3;
00040
00041 zero=Mp->zero;
00042 }
00043
00044
00045
00046 barelq3d::~barelq3d (void)
00047 {
00048 long i;
00049
00050 for (i=0;i<nb;i++){
00051 delete [] nip[i];
00052 delete [] intordsm[i];
00053 }
00054 delete [] nip;
00055 delete [] intordsm;
00056
00057 delete [] cncomp;
00058 delete [] ncomp;
00059 }
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 void barelq3d::dirvect (vector &s,vector &x,vector &y,vector &z)
00075 {
00076 double d;
00077
00078 s[0]=x[1]-x[0];
00079 s[1]=y[1]-y[0];
00080 s[2]=z[1]-z[0];
00081
00082 d=sqrt(s[0]*s[0]+s[1]*s[1]+s[2]*s[2]);
00083
00084 if (d<zero)
00085 print_err("zero norm of direction vector",__FILE__,__LINE__,__func__);
00086
00087 s[0]/=d;
00088 s[1]/=d;
00089 s[2]/=d;
00090 }
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 double barelq3d::approx (double xi,vector &nodval)
00105 {
00106 double f;
00107 vector bf(nne);
00108
00109 bf_quad_1d (bf.a,xi);
00110
00111 scprd (bf,nodval,f);
00112
00113 return f;
00114 }
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 void barelq3d::geom_matrix (matrix &gm, vector &x,vector &s,double xi,double &jac)
00133 {
00134 vector dx(nne);
00135
00136 dx_bf_quad_1d (dx.a,xi);
00137 derivatives_1d (dx, jac, x, xi);
00138
00139 fillm (0.0, gm);
00140
00141 gm[0][0]= dx[0]*s[0];
00142 gm[0][1]= dx[0]*s[1];
00143 gm[0][2]= dx[0]*s[2];
00144 gm[0][3]= dx[1]*s[0];
00145 gm[0][4]= dx[1]*s[1];
00146 gm[0][5]= dx[1]*s[2];
00147 gm[0][6]= dx[2]*s[0];
00148 gm[0][7]= dx[2]*s[1];
00149 gm[0][8]= dx[2]*s[2];
00150 }
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 void barelq3d::giveloccoord(vector &x,vector &y,vector &z,vector &lx)
00166 {
00167 double l;
00168
00169 lx[0] = 0.0;
00170 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]));
00171 lx[1] = l;
00172 l = sqrt((x[2]-x[0])*(x[2]-x[0]) + (y[2]-y[0])*(y[2]-y[0]) + (z[2]-z[0])*(z[2]-z[0]));
00173 lx[2] = l;
00174 }
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 void barelq3d::give_glob_loc_tmat(vector &x, vector &y, matrix &tmat)
00193 {
00194 double c, s, l;
00195
00196 fillm(0.0, tmat);
00197 l = sqrt(sqr(x[1]-x[0]) + sqr(y[1]-y[0]));
00198 c = (x[1] - x[0]) / l;
00199 s = (y[1] - y[0]) / l;
00200 tmat[0][0] = tmat[2][2] = tmat[4][4] = c;
00201 tmat[0][1] = tmat[2][3] = tmat[4][5] = s;
00202 tmat[1][0] = tmat[3][2] = tmat[5][4] = -s;
00203 tmat[1][1] = tmat[3][3] = tmat[5][5] = c;
00204 }
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222 void barelq3d::give_loc_glob_tmat(vector &x, vector &y, matrix &tmat)
00223 {
00224 double c, s, l;
00225
00226 fillm(0.0, tmat);
00227 l = sqrt(sqr(x[1]-x[0]) + sqr(y[1]-y[0]));
00228 c = (x[1] - x[0]) / l;
00229 s = (y[1] - y[0]) / l;
00230 tmat[0][0] = tmat[2][2] = tmat[4][4] = c;
00231 tmat[0][1] = tmat[2][3] = tmat[4][5] = -s;
00232 tmat[1][0] = tmat[3][2] = tmat[5][4] = s;
00233 tmat[1][1] = tmat[3][3] = tmat[5][5] = c;
00234 }
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 void barelq3d::transf_matrix (ivector &nodes,matrix &tmat)
00248 {
00249 long i;
00250 fillm (0.0,tmat);
00251
00252 for (i=0;i<tmat.m;i++){
00253 tmat[i][i]=1.0;
00254 }
00255
00256
00257 for (i=0;i<nodes.n;i++){
00258 if (Mt->nodes[nodes[i]].transf>0){
00259 tmat[i*3+0][i*3]=Mt->nodes[nodes[i]].e1[0];
00260 tmat[i*3+1][i*3]=Mt->nodes[nodes[i]].e1[1];
00261 tmat[i*3+2][i*3]=Mt->nodes[nodes[i]].e1[2];
00262
00263 tmat[i*3+0][i*3+1]=Mt->nodes[nodes[i]].e2[0];
00264 tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e2[1];
00265 tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e2[2];
00266
00267 tmat[i*3+0][i*3+2]=Mt->nodes[nodes[i]].e3[0];
00268 tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e3[1];
00269 tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e3[2];
00270 }
00271 }
00272 }
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 void barelq3d::stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &lx,vector &s)
00290 {
00291 long i,ipp;
00292 double a,xi,jac;
00293 vector w,gp,t(nne);
00294 matrix gm(tncomp,ndofe),d(tncomp,tncomp);
00295
00296
00297 fillm (0.0,sm);
00298
00299 Mc->give_area (eid,a);
00300
00301 allocv (intordsm[0][0],w);
00302 allocv (intordsm[0][0],gp);
00303 gauss_points (gp.a,w.a,intordsm[0][0]);
00304
00305
00306 ipp=Mt->elements[eid].ipp[ri][ci];
00307
00308 for (i=0;i<intordsm[0][0];i++){
00309 xi=gp[i];
00310
00311
00312 geom_matrix (gm,lx,s,xi,jac);
00313
00314 Mm->matstiff (d,ipp);
00315 jac*=a*w[i];
00316
00317 bdbjac (sm,gm,d,gm,jac);
00318
00319 ipp++;
00320 }
00321
00322 }
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337 void barelq3d::res_stiffness_matrix (long eid,matrix &sm)
00338 {
00339 long transf;
00340 ivector nodes(nne);
00341 vector gx(nne),gy(nne),gz(nne),x(nne),s(nne);
00342
00343
00344 Mt->give_node_coord3d (gx,gy,gz,eid);
00345
00346
00347 giveloccoord (gx,gy,gz,x);
00348
00349 dirvect (s,gx,gy,gz);
00350
00351
00352 stiffness_matrix (eid,0,0,sm,x,s);
00353
00354
00355
00356 Mt->give_elemnodes (eid,nodes);
00357 transf = Mt->locsystems (nodes);
00358 if (transf>0){
00359 matrix tmat (ndofe,ndofe);
00360 transf_matrix (nodes,tmat);
00361 glmatrixtransf (sm,tmat);
00362 }
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 void barelq3d::mass_matrix (long eid,matrix &mm)
00377 {
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414 }
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426 void barelq3d::res_mainip_strains (long lcid,long eid)
00427 {
00428 vector aux,gx(nne),gy(nne),gz(nne),x(nne),s(3),r(ndofe);
00429 ivector nodes(nne);
00430 matrix tmat;
00431
00432
00433 Mt->give_node_coord3d (gx,gy,gz,eid);
00434
00435 giveloccoord (gx,gy,gz,x);
00436
00437 dirvect (s,gx,gy,gz);
00438
00439 Mt->give_elemnodes (eid,nodes);
00440
00441 eldispl (lcid,eid,r.a);
00442
00443
00444 long transf = Mt->locsystems (nodes);
00445 if (transf>0){
00446 allocv (ndofe,aux);
00447 allocm (ndofe,ndofe,tmat);
00448 transf_matrix (nodes,tmat);
00449
00450 lgvectortransf (aux,r,tmat);
00451 copyv (aux,r);
00452 destrv (aux);
00453 destrm (tmat);
00454 }
00455
00456
00457 mainip_strains (lcid,eid,0,0,x,s,r);
00458 }
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475 void barelq3d::mainip_strains (long lcid,long eid,long ri,long ci,vector &lx,vector &s,vector &r)
00476 {
00477 long i,ii,ipp;
00478 double xi,jac;
00479 vector gp,w,eps;
00480 matrix gm;
00481
00482 for (ii=0;ii<nb;ii++){
00483 if (intordsm[ii][ii]==0) continue;
00484
00485 allocv (intordsm[ii][ii],gp);
00486 allocv (intordsm[ii][ii],w);
00487 allocv (ncomp[ii],eps);
00488 allocm (ncomp[ii],ndofe,gm);
00489
00490 gauss_points (gp.a,w.a,intordsm[ii][ii]);
00491
00492
00493 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00494
00495 for (i=0;i<intordsm[ii][ii];i++){
00496 xi=gp[i];
00497
00498
00499 geom_matrix (gm,lx,s,xi,jac);
00500 mxv (gm,r,eps);
00501
00502 Mm->storestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00503 ipp++;
00504 }
00505
00506 destrm (gm); destrv (eps); destrv (w); destrv (gp);
00507 }
00508 }
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519 void barelq3d::nod_strains_ip (long lcid,long eid,long ri,long ci)
00520 {
00521 long i,j;
00522 ivector ipnum(nne),nod(nne);
00523 vector eps(tncomp);
00524
00525
00526 nodipnum (eid,ri,ci,ipnum);
00527
00528
00529 Mt->give_elemnodes (eid,nod);
00530
00531 for (i=0;i<nne;i++){
00532
00533 Mm->givestrain (lcid,ipnum[i],eps);
00534
00535
00536 j=nod[i];
00537 Mt->nodes[j].storestrain (lcid,0,eps);
00538 }
00539
00540 }
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555 void barelq3d::nod_strains_comp (long lcid,long eid,double **stra)
00556 {
00557 long i,j;
00558 double jac;
00559 ivector nodes(nne);
00560 vector gx(nne),gy(nne),gz(nne),x(nne),s(3),nxi(nne),r(ndofe),eps(tncomp),aux;
00561 matrix tmat,gm(tncomp,ndofe);
00562
00563
00564 Mt->give_node_coord3d (gx,gy,gz,eid);
00565
00566 giveloccoord (gx,gy,gz,x);
00567
00568 dirvect (s,gx,gy,gz);
00569
00570 Mt->give_elemnodes (eid,nodes);
00571
00572 eldispl (lcid,eid,r.a);
00573
00574
00575 long transf = Mt->locsystems (nodes);
00576 if (transf>0){
00577 allocv (ndofe,aux);
00578 allocm (ndofe,ndofe,tmat);
00579 transf_matrix (nodes,tmat);
00580
00581 lgvectortransf (aux,r,tmat);
00582 copyv (aux,r);
00583 destrv (aux);
00584 destrm (tmat);
00585 }
00586
00587
00588 nodecoord (nxi);
00589
00590
00591 for (i=0;i<nne;i++){
00592
00593 geom_matrix (gm,x,s,nxi[i],jac);
00594
00595 mxv (gm,r,eps);
00596
00597 for (j=0;j<eps.n;j++){
00598 stra[i][j]=eps[j];
00599 }
00600 }
00601 }
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 void barelq3d::res_allip_strains (long lcid,long eid)
00614 {
00615
00616 allip_strains (lcid,eid,0,0);
00617 }
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630 void barelq3d::allip_strains (long lcid,long eid,long ri,long ci)
00631 {
00632
00633 res_mainip_strains (lcid,eid);
00634 }
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 void barelq3d::strains (long lcid,long eid,long ri,long ci)
00648 {
00649
00650
00651 vector coord,eps;
00652
00653 switch (Mm->stra.tape[eid]){
00654 case nowhere:{
00655 break;
00656 }
00657 case intpts:{
00658 allip_strains (lcid,eid,ri,ci);
00659 break;
00660 }
00661 case enodes:{
00662 nod_strains_ip (lcid,eid,ri,ci);
00663 break;
00664 }
00665 case userdefined:{
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691 break;
00692 }
00693 default:
00694 print_err("unknown strain point is required",__FILE__,__LINE__,__func__);
00695 }
00696
00697 }
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710 void barelq3d::nodecoord (vector &xi)
00711 {
00712 xi[0] = -1.0;
00713 xi[1] = 1.0;
00714 xi[2] = 0.0;
00715 }
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730 void barelq3d::nodipnum (long eid,long ri,long ci,ivector &ipnum)
00731 {
00732 long i;
00733
00734 i=Mt->elements[eid].ipp[ri][ci];
00735
00736 ipnum[0]=i;
00737 ipnum[1]=i+2;
00738 ipnum[2]=i+1;
00739 }
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751 void barelq3d::res_mainip_stresses (long lcid,long eid)
00752 {
00753 mainip_stresses (lcid,eid,0,0);
00754 }
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768 void barelq3d::mainip_stresses (long lcid,long eid,long ri,long ci)
00769 {
00770 long i,ii,jj,ipp;
00771 double xi;
00772 vector gp,w,eps,epst,epstt,sig,auxsig;
00773 matrix d(tncomp,tncomp);
00774
00775 for (ii=0;ii<nb;ii++){
00776
00777 allocv (ncomp[ii],sig);
00778 allocv (ncomp[ii],auxsig);
00779 allocv (intordsm[ii][ii],gp);
00780 allocv (intordsm[ii][ii],w);
00781
00782 gauss_points (gp.a,w.a,intordsm[ii][ii]);
00783 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00784
00785 for (i=0;i<intordsm[ii][ii];i++){
00786 xi=gp[i];
00787
00788 Mm->matstiff (d,ipp);
00789
00790 fillv (0.0,sig);
00791 for (jj=0;jj<nb;jj++){
00792 allocv (ncomp[jj],eps);
00793
00794
00795 Mm->givestrain (lcid,ipp,cncomp[jj],ncomp[jj],eps);
00796
00797
00798 mxv (d,eps,auxsig);
00799
00800 addv (auxsig,sig,sig);
00801
00802 destrv (eps);
00803 }
00804
00805
00806 Mm->storestress (lcid,ipp,cncomp[ii],ncomp[ii],sig);
00807
00808 ipp++;
00809 }
00810
00811 destrv (w); destrv (gp); destrv (auxsig); destrv (sig);
00812 }
00813 }
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 void barelq3d::nod_stresses_ip (long lcid,long eid,long ri,long ci)
00827 {
00828 long i,j;
00829 ivector ipnum(nne),nod(nne);
00830 vector sig(tncomp);
00831
00832
00833 nodipnum (eid,ri,ci,ipnum);
00834
00835
00836 Mt->give_elemnodes (eid,nod);
00837
00838 for (i=0;i<nne;i++){
00839
00840 Mm->givestress (lcid,ipnum[i],sig);
00841
00842
00843 j=nod[i];
00844 Mt->nodes[j].storestress (lcid,0,sig);
00845 }
00846
00847 }
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864 void barelq3d::nod_stresses_comp (long lcid,long eid,long ri,long ci,double **stra,double **stre)
00865 {
00866 long i,j,ipp;
00867 vector eps(tncomp),sig(tncomp);
00868 matrix d(tncomp,tncomp);
00869
00870
00871 ipp=Mt->elements[eid].ipp[ri][ci];
00872
00873 Mm->matstiff (d,ipp);
00874
00875
00876 for (i=0;i<nne;i++){
00877 for (j=0;j<eps.n;j++){
00878 eps[j]=stra[i][j];
00879 }
00880 mxv (d,eps,sig);
00881 for (j=0;j<eps.n;j++){
00882 stre[i][j]=sig[j];
00883 }
00884 }
00885 }
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897 void barelq3d::res_allip_stresses (long lcid,long eid)
00898 {
00899
00900 allip_stresses (lcid,eid,0,0);
00901 }
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914 void barelq3d::allip_stresses (long lcid,long eid,long ri,long ci)
00915 {
00916
00917 res_mainip_stresses (lcid,eid);
00918 }
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930 void barelq3d::stresses (long lcid,long eid,long ri,long ci)
00931 {
00932 vector coord,sig;
00933
00934 switch (Mm->stre.tape[eid]){
00935 case nowhere:{
00936 break;
00937 }
00938 case intpts:{
00939 allip_stresses (lcid,eid,ri,ci);
00940 break;
00941 }
00942 case enodes:{
00943 nod_stresses_ip (lcid,eid,ri,ci);
00944 break;
00945 }
00946 case userdefined:{
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967 break;
00968 }
00969 default:
00970 print_err("unknown stress point is required",__FILE__,__LINE__,__func__);
00971 }
00972
00973 }
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986 void barelq3d::nod_eqother_ip (long lcid,long eid,long ri,long ci)
00987 {
00988 long i,j,ncompo;
00989 ivector ipnum(nne),nod(nne);
00990 vector eqother;
00991
00992
00993 nodipnum (eid,ri,ci,ipnum);
00994
00995
00996 Mt->give_elemnodes (eid,nod);
00997
00998 for (i=0;i<nne;i++){
00999
01000
01001
01002 ncompo = Mm->givencompeqother (ipnum[i],0);
01003 allocv (ncompo,eqother);
01004 Mm->giveeqother (ipnum[i],0,ncompo,eqother.a);
01005
01006
01007 j=nod[i];
01008 Mt->nodes[j].storeother (lcid,0,ncompo,eqother);
01009
01010 destrv (eqother);
01011 }
01012 }
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027 void barelq3d::load_matrix (long eid,matrix &lm)
01028 {
01029
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
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075 void barelq3d::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &s)
01076 {
01077 integratedquant iq;
01078 iq=locstress;
01079
01080
01081 compute_nlstress (lcid,eid,ri,ci);
01082
01083
01084 elem_integration (iq,lcid,eid,ri,ci,ifor,x,s);
01085 }
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103 void barelq3d::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &s)
01104 {
01105 integratedquant iq;
01106 iq=nonlocstress;
01107
01108
01109 compute_nonloc_nlstress (lcid,eid,ri,ci);
01110
01111
01112 elem_integration (iq,lcid,eid,ri,ci,ifor,x,s);
01113 }
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131 void barelq3d::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &s)
01132 {
01133 integratedquant iq;
01134 iq=stressincr;
01135
01136
01137 compute_nlstressincr (lcid,eid,ri,ci);
01138
01139
01140 elem_integration (iq,lcid,eid,ri,ci,ifor,x,s);
01141 }
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158 void barelq3d::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &s)
01159 {
01160 integratedquant iq;
01161 iq=eigstress;
01162
01163
01164 if (Mp->eigstrcomp)
01165 compute_eigstress (lcid,eid,ri,ci);
01166
01167
01168 elem_integration (iq,lcid,eid,ri,ci,nfor,x,s);
01169 }
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185 void barelq3d::res_internal_forces (long lcid,long eid,vector &ifor)
01186 {
01187 vector gx(nne),gy(nne),gz(nne),x(nne),s(3);
01188 long transf;
01189 ivector nodes (nne);
01190
01191
01192 Mt->give_node_coord3d (gx,gy,gz,eid);
01193
01194 giveloccoord (gx,gy,gz,x);
01195
01196 dirvect (s,gx,gy,gz);
01197
01198 internal_forces (lcid,eid,0,0,ifor,x,s);
01199
01200
01201
01202 Mt->give_elemnodes (eid,nodes);
01203 transf = Mt->locsystems (nodes);
01204 if (transf>0){
01205 matrix tmat (ndofe,ndofe);
01206 vector v(ndofe);
01207 transf_matrix (nodes,tmat);
01208 glvectortransf (ifor,v,tmat);
01209 copyv (v,ifor);
01210 }
01211 }
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226 void barelq3d::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
01227 {
01228 vector gx(nne),gy(nne),gz(nne),x(nne),s(3);
01229 long transf;
01230 ivector nodes(nne);
01231
01232
01233 Mt->give_node_coord3d (gx,gy,gz,eid);
01234
01235 giveloccoord (gx,gy,gz,x);
01236
01237 dirvect (s,gx,gy,gz);
01238
01239 nonloc_internal_forces (lcid,eid,0,0,ifor,x,s);
01240
01241
01242
01243 Mt->give_elemnodes (eid,nodes);
01244 transf = Mt->locsystems (nodes);
01245 if (transf>0){
01246 matrix tmat (ndofe,ndofe);
01247 vector v(ndofe);
01248 transf_matrix (nodes,tmat);
01249 glvectortransf (ifor,v,tmat);
01250 copyv (v,ifor);
01251 }
01252 }
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267 void barelq3d::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01268 {
01269 vector gx(nne),gy(nne),gz(nne),x(nne),s(3);
01270 long transf;
01271 ivector nodes(nne);
01272
01273
01274 Mt->give_node_coord3d (gx,gy,gz,eid);
01275
01276 giveloccoord (gx,gy,gz,x);
01277
01278 dirvect (s,gx,gy,gz);
01279
01280 incr_internal_forces (lcid,eid,0,0,ifor,x,s);
01281
01282
01283
01284 Mt->give_elemnodes (eid,nodes);
01285 transf = Mt->locsystems (nodes);
01286 if (transf>0){
01287 matrix tmat (ndofe,ndofe);
01288 vector v(ndofe);
01289 transf_matrix (nodes,tmat);
01290 glvectortransf (ifor,v,tmat);
01291 copyv (v,ifor);
01292 }
01293 }
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308 void barelq3d::res_eigstrain_forces (long lcid,long eid,vector &ifor)
01309 {
01310 vector gx(nne),gy(nne),gz(nne),x(nne),s(3);
01311 long transf;
01312 ivector nodes(nne);
01313
01314
01315 Mt->give_node_coord3d (gx,gy,gz,eid);
01316
01317 giveloccoord (gx,gy,gz,x);
01318
01319 dirvect (s,gx,gy,gz);
01320
01321 eigstrain_forces (lcid,eid,0,0,ifor,x,s);
01322
01323
01324
01325 Mt->give_elemnodes (eid,nodes);
01326 transf = Mt->locsystems (nodes);
01327 if (transf>0){
01328 matrix tmat (ndofe,ndofe);
01329 vector v(ndofe);
01330 transf_matrix (nodes,tmat);
01331 glvectortransf (ifor,v,tmat);
01332 copyv (v,ifor);
01333 }
01334 }
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353 void barelq3d::elem_integration (integratedquant iq, long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &s)
01354 {
01355 long i,ii,ipp;
01356 double xi,jac,area;
01357 vector w,gp,ipv,contr(ndofe);
01358 matrix gm;
01359
01360
01361 Mc->give_area (eid,area);
01362 fillv (0.0,nv);
01363
01364 for (ii=0;ii<nb;ii++){
01365 allocv (intordsm[ii][ii],gp);
01366 allocv (intordsm[ii][ii],w);
01367 allocm (ncomp[ii],ndofe,gm);
01368 allocv (ncomp[ii],ipv);
01369
01370 gauss_points (gp.a,w.a,intordsm[ii][ii]);
01371 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01372
01373 for (i=0;i<intordsm[ii][ii];i++){
01374 xi=gp[i];
01375
01376
01377 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
01378
01379
01380 geom_matrix (gm,x,s,xi,jac);
01381
01382
01383 mtxv (gm,ipv,contr);
01384
01385 cmulv (jac*w[i]*area,contr);
01386
01387
01388 addv(contr, nv, nv);
01389
01390 ipp++;
01391 }
01392 }
01393 destrv (ipv); destrm (gm); destrv (w); destrv (gp);
01394
01395 }
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408 void barelq3d::compute_nlstress (long lcid,long eid,long ri,long ci)
01409 {
01410 long i,ii,ipp;
01411
01412 for (ii=0;ii<nb;ii++){
01413 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01414
01415 for (i=0;i<intordsm[ii][ii];i++){
01416
01417 if (Mp->strcomp==1)
01418 Mm->computenlstresses (ipp);
01419
01420 ipp++;
01421 }
01422 }
01423 }
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436 void barelq3d::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01437 {
01438 long i,ii,ipp;
01439
01440 for (ii=0;ii<nb;ii++){
01441 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01442
01443 for (i=0;i<intordsm[ii][ii];i++){
01444
01445 if (Mp->strcomp==1)
01446 Mm->computenlstressesincr (ipp);
01447
01448 ipp++;
01449 }
01450 }
01451 }
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465 void barelq3d::local_values (long lcid,long eid,long ri,long ci)
01466 {
01467 long i,ii,ipp;
01468
01469 for (ii=0;ii<nb;ii++){
01470 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01471
01472 for (i=0;i<intordsm[ii][ii];i++){
01473
01474 if (Mp->strcomp==1)
01475 Mm->computenlstresses (ipp);
01476
01477 ipp++;
01478 }
01479 }
01480 }
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493 void barelq3d::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01494 {
01495 long i,ii,ipp;
01496
01497 for (ii=0;ii<nb;ii++){
01498 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01499
01500 for (i=0;i<intordsm[ii][ii];i++){
01501
01502 if (Mp->strcomp==1)
01503 Mm->compnonloc_nlstresses (ipp);
01504
01505 ipp++;
01506 }
01507 }
01508 }
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521 void barelq3d::compute_eigstress (long lcid,long eid,long ri,long ci)
01522 {
01523 long i,ipp;
01524 vector eigstr(tncomp),sig(tncomp);
01525 matrix d(tncomp,tncomp);
01526
01527 ipp=Mt->elements[eid].ipp[ri][ci];
01528
01529 for (i=0;i<intordsm[0][0];i++){
01530
01531 Mm->giveeigstrain (ipp,eigstr);
01532
01533
01534 Mm->matstiff (d,ipp);
01535
01536 mxv (d,eigstr,sig);
01537
01538 Mm->storeeigstress (ipp,sig);
01539
01540 ipp++;
01541 }
01542 }
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555 void barelq3d::res_mass_matrix (long eid,matrix &mm)
01556 {
01557 vector gx(nne),gy(nne),gz(nne),x(nne);
01558
01559
01560 Mt->give_node_coord2d (gx,gy,eid);
01561 giveloccoord (gx,gy,gz,x);
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576 }
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602 void barelq3d::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
01603 {
01604 long i, j, k, ipp;
01605 long ii, jj, nv = nodval.n;
01606 double xi, ipval;
01607 vector w, gp, anv(3);
01608 long nstra, nstre, ncompstr, ncompeqother;
01609 long idstra, idstre, idoth, idic;
01610 inictype ict;
01611
01612 nstra = idstra = nstre = idstre = idoth = idic = 0;
01613
01614 ict = ictn[0];
01615 for (i=0; i<nne; i++)
01616 {
01617 if (ictn[i] != ict)
01618 {
01619 print_err("Incompatible types of initial conditions on element %ld\n"
01620 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
01621 abort();
01622 }
01623 }
01624 for (j = 0; j < nv; j++)
01625 {
01626 for(i = 0; i < nne; i++)
01627 anv[i] = nodval[i][j];
01628 for (ii = 0; ii < nb; ii++)
01629 {
01630 for (jj = 0; jj < nb; jj++)
01631 {
01632 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01633 if (intordsm[ii][jj] == 0)
01634 continue;
01635 allocv (intordsm[ii][jj],gp);
01636 allocv (intordsm[ii][jj],w);
01637 gauss_points (gp.a,w.a,intordsm[ii][jj]);
01638 for (k = 0; k < intordsm[ii][jj]; k++)
01639 {
01640 xi=gp[k];
01641
01642 ipval = approx (xi,anv);
01643 ncompstr = Mm->ip[ipp].ncompstr;
01644 ncompeqother = Mm->ip[ipp].ncompeqother;
01645 if ((ictn[0] & inistrain) && (j < ncompstr))
01646 {
01647 Mm->ip[ipp].strain[idstra] += ipval;
01648 ipp++;
01649 continue;
01650 }
01651 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01652 {
01653 Mm->ip[ipp].stress[idstre] += ipval;
01654 ipp++;
01655 continue;
01656 }
01657 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
01658 {
01659 Mm->ip[ipp].eqother[idoth] += ipval;
01660 ipp++;
01661 continue;
01662 }
01663 if ((ictn[0] & inicond) && (j < nv))
01664 {
01665 if (Mm->ic[ipp] == NULL)
01666 {
01667 Mm->ic[ipp] = new double[nv-j];
01668 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
01669 }
01670 Mm->ic[ipp][idic] += ipval;
01671 ipp++;
01672 continue;
01673 }
01674 ipp++;
01675 }
01676 destrv (gp); destrv (w);
01677 }
01678 }
01679 ipp=Mt->elements[eid].ipp[ri][ci];
01680 ncompstr = Mm->ip[ipp].ncompstr;
01681 ncompeqother = Mm->ip[ipp].ncompeqother;
01682 if ((ictn[0] & inistrain) && (j < ncompstr))
01683 {
01684 nstra++;
01685 idstra++;
01686 continue;
01687 }
01688 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01689 {
01690 nstre++;
01691 idstre++;
01692 continue;
01693 }
01694 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
01695 {
01696 idoth++;
01697 continue;
01698 }
01699 if ((ictn[0] & inicond) && (j < nv))
01700 {
01701 idic++;
01702 continue;
01703 }
01704 }
01705 }
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719 void barelq3d::intpointval (long eid,vector &nodval,vector &ipval)
01720 {
01721 long ii,jj,i,k,ipp;
01722 double xi;
01723 vector w,gp;
01724
01725 k=0;
01726 for (ii = 0; ii < nb; ii++)
01727 {
01728 for (jj = 0; jj < nb; jj++)
01729 {
01730 ipp=Mt->elements[eid].ipp[ii][jj];
01731 if (intordsm[ii][jj] == 0)
01732 continue;
01733 allocv (intordsm[ii][jj],gp);
01734 allocv (intordsm[ii][jj],w);
01735 gauss_points (gp.a,w.a,intordsm[ii][jj]);
01736 for (i = 0; i < intordsm[ii][jj]; i++)
01737 {
01738 xi=gp[i];
01739
01740 ipval[k] = approx (xi,nodval);
01741 k++;
01742 ipp++;
01743 }
01744 destrv (gp); destrv (w);
01745 }
01746 }
01747 }
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761 void barelq3d::intpointval2 (long eid,vector &nodval,vector &ipval)
01762 {
01763 long ii,jj,i,k;
01764 double xi;
01765 vector w,gp;
01766 vector modnodval(Bar2d->nne);
01767
01768 for (i=0;i<Bar2d->nne;i++){
01769 modnodval[i]=nodval[i];
01770 }
01771
01772 k=0;
01773 for (ii = 0; ii < nb; ii++){
01774 for (jj = 0; jj < nb; jj++){
01775 if (intordsm[ii][jj] == 0)
01776 continue;
01777
01778 allocv (intordsm[ii][jj],gp);
01779 allocv (intordsm[ii][jj],w);
01780
01781 gauss_points (gp.a,w.a,intordsm[ii][jj]);
01782
01783 for (i = 0; i < intordsm[ii][jj]; i++){
01784 xi=gp[i];
01785
01786
01787 ipval[k] = Bar2d->approx (xi,modnodval);
01788
01789 k++;
01790 }
01791 destrv (gp); destrv (w);
01792 }
01793 }
01794 }
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811 void barelq3d::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
01812 {
01813 long ii,jj,i,aipp;
01814 double xi;
01815 vector w,gp;
01816 vector x(nne),y(nne),z(nne);
01817
01818 Mt->give_node_coord3d (x,y,z,eid);
01819 for (ii = 0; ii < nb; ii++){
01820 for (jj = 0; jj < nb; jj++){
01821 if (intordsm[ii][jj] == 0)
01822 continue;
01823
01824 aipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01825 allocv (intordsm[ii][jj],gp);
01826 allocv (intordsm[ii][jj],w);
01827 gauss_points (gp.a,w.a,intordsm[ii][jj]);
01828
01829 for (i=0;i<intordsm[ii][jj];i++){
01830 xi=gp[i];
01831 if (aipp==ipp){
01832 coord[0]=approx (xi,x);
01833 coord[1]=approx (xi,y);
01834 coord[2]=approx (xi,z);
01835 }
01836 aipp++;
01837 }
01838 destrv (gp); destrv (w);
01839 }
01840 }
01841 }