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