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