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