00001 #include "plelemlq.h"
00002 #include "global.h"
00003 #include "globmat.h"
00004 #include "genfile.h"
00005 #include "adaptivity.h"
00006 #include "node.h"
00007 #include "element.h"
00008 #include "intpoints.h"
00009 #include "loadcase.h"
00010 #include "gadaptivity.h"
00011 #include <stdlib.h>
00012 #include <math.h>
00013
00014
00015
00016 planeelemlq::planeelemlq (void)
00017 {
00018 long i,j;
00019
00020
00021 nne=4;
00022
00023 ndofe=8;
00024
00025 tncomp=3;
00026
00027 gncomp=4;
00028
00029 napfun=2;
00030
00031 intordmm=2;
00032
00033 ned=4;
00034
00035 nned=2;
00036
00037 intordb=2;
00038
00039
00040 nb=2;
00041
00042
00043 ncomp = new long [nb];
00044 ncomp[0]=2;
00045 ncomp[1]=1;
00046
00047
00048 cncomp = new long [nb];
00049 cncomp[0]=0;
00050 cncomp[1]=2;
00051
00052
00053
00054 nip = new long* [nb];
00055 intordsm = new long* [nb];
00056 for (i=0;i<nb;i++){
00057 nip[i] = new long [nb];
00058 intordsm[i] = new long [nb];
00059 }
00060
00061 nip[0][0]=4; nip[0][1]=0;
00062 nip[1][0]=0; nip[1][1]=1;
00063
00064
00065 tnip=0;
00066 for (i=0;i<nb;i++){
00067 for (j=0;j<nb;j++){
00068 tnip+=nip[i][j];
00069 }
00070 }
00071
00072 intordsm[0][0]=2; intordsm[0][1]=0;
00073 intordsm[1][0]=0; intordsm[1][1]=1;
00074
00075 }
00076
00077 planeelemlq::~planeelemlq (void)
00078 {
00079 long i;
00080
00081 for (i=0;i<nb;i++){
00082 delete [] nip[i];
00083 delete [] intordsm[i];
00084 }
00085 delete [] nip;
00086 delete [] intordsm;
00087
00088 delete [] cncomp;
00089 delete [] ncomp;
00090 }
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101 double planeelemlq::approx (double xi,double eta,vector &nodval)
00102 {
00103 double f;
00104 vector bf(nne);
00105
00106 bf_lin_4_2d (bf.a,xi,eta);
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 void planeelemlq::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
00126 {
00127 long i,j,ii;
00128 double xi,eta;
00129 vector x(nne),y(nne),w(intordsm[ri][ci]),gp(intordsm[ri][ci]);
00130
00131 gauss_points (gp.a,w.a,intordsm[ri][ci]);
00132 Mt->give_node_coord2d (x,y,eid);
00133 ii=Mt->elements[eid].ipp[ri][ci];
00134
00135 for (i=0;i<intordsm[ri][ci];i++){
00136 xi=gp[i];
00137 for (j=0;j<intordsm[ri][ci];j++){
00138 eta=gp[j];
00139 if (ii==ipp){
00140 coord[0]=approx (xi,eta,x);
00141 coord[1]=approx (xi,eta,y);
00142 coord[2]=0.0;
00143 }
00144 ii++;
00145 }
00146 }
00147 }
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 void planeelemlq::ipcoordblock (long eid,long ri,long ci,double **coord)
00160 {
00161 long i,j,k;
00162 double xi,eta;
00163 vector x(nne),y(nne),w(intordsm[ri][ci]),gp(intordsm[ri][ci]);
00164
00165 gauss_points (gp.a,w.a,intordsm[ri][ci]);
00166 Mt->give_node_coord2d (x,y,eid);
00167
00168 k=0;
00169 for (i=0;i<intordsm[ri][ci];i++){
00170 xi=gp[i];
00171 for (j=0;j<intordsm[ri][ci];j++){
00172 eta=gp[j];
00173
00174 coord[k][0]=approx (xi,eta,x);
00175 coord[k][1]=approx (xi,eta,y);
00176 coord[k][2]=0.0;
00177
00178 k++;
00179 }
00180 }
00181 }
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 void planeelemlq::bf_matrix (matrix &n,double xi,double eta)
00192 {
00193 long i,j,k;
00194 vector bf(nne);
00195
00196 bf_lin_4_2d (bf.a,xi,eta);
00197
00198 fillm (0.0,n);
00199
00200 j=0; k=1;
00201 for (i=0;i<nne;i++){
00202 n[0][j]=bf[i];
00203 n[1][k]=bf[i];
00204 j+=2; k+=2;
00205 }
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220 void planeelemlq::geom_matrix (matrix &gm,vector &x,vector &y,double xi,double eta,double &jac)
00221 {
00222 long i,i1,i2;
00223 vector dx(nne),dy(nne);
00224
00225 dx_bf_lin_4_2d (dx.a,eta);
00226 dy_bf_lin_4_2d (dy.a,xi);
00227
00228 derivatives_2d (dx,dy,jac,x,y,xi,eta);
00229
00230 fillm (0.0,gm);
00231
00232 i1=0; i2=1;
00233 for (i=0;i<nne;i++){
00234 gm[0][i1]=dx[i];
00235 gm[1][i2]=dy[i];
00236 gm[2][i1]=dy[i];
00237 gm[2][i2]=dx[i];
00238 i1+=2; i2+=2;
00239 }
00240
00241 }
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 void planeelemlq::geom_matrix_block (matrix &gm,long ri,vector &x,vector &y,double xi,double eta,double &jac)
00255 {
00256 long i,i1,i2;
00257 vector dx(nne),dy(nne);
00258
00259 dx_bf_lin_4_2d (dx.a,eta);
00260 dy_bf_lin_4_2d (dy.a,xi);
00261
00262 derivatives_2d (dx,dy,jac,x,y,xi,eta);
00263
00264 fillm (0.0,gm);
00265
00266 if (ri==0){
00267 i1=0; i2=1;
00268 for (i=0;i<nne;i++){
00269 gm[0][i1]=dx[i];
00270 gm[1][i2]=dy[i];
00271 i1+=2; i2+=2;
00272 }
00273 }
00274 if (ri==1){
00275 i1=0; i2=1;
00276 for (i=0;i<nne;i++){
00277 gm[0][i1]=dy[i];
00278 gm[0][i2]=dx[i];
00279 i1+=2; i2+=2;
00280 }
00281 }
00282 }
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 void planeelemlq::bvectors (vector &x,vector &y,double xi,double eta,double &jac,
00296 vector &b11,vector &b12,vector &b21,vector &b22)
00297 {
00298 vector dx(nne),dy(nne);
00299
00300 dx_bf_lin_4_2d (dx.a,eta);
00301 dy_bf_lin_4_2d (dy.a,xi);
00302
00303 derivatives_2d (dx,dy,jac,x,y,xi,eta);
00304
00305 fillv (0.0,b11);
00306 fillv (0.0,b12);
00307 fillv (0.0,b21);
00308 fillv (0.0,b22);
00309
00310
00311 b11[0]=dx[0]; b11[2]=dx[1]; b11[4]=dx[2]; b11[6]=dx[3];
00312
00313 b12[0]=dy[0]; b12[2]=dy[1]; b12[4]=dy[2]; b12[6]=dy[3];
00314
00315 b21[1]=dx[0]; b21[3]=dx[1]; b21[5]=dx[2]; b21[7]=dx[3];
00316
00317 b22[1]=dy[0]; b22[3]=dy[1]; b22[5]=dy[2]; b22[7]=dy[3];
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331 void planeelemlq::gngeom_matrix (matrix &gm,vector &r,vector &x,vector &y,double xi,double eta,double &jac)
00332 {
00333 long i;
00334 vector b11(ndofe),b12(ndofe),b21(ndofe),b22(ndofe),av(ndofe);
00335 matrix am(ndofe,ndofe);
00336
00337 fillm (0.0,gm);
00338
00339 bvectors (x,y,xi,eta,jac,b11,b12,b21,b22);
00340
00341
00342
00343
00344
00345
00346 for (i=0;i<ndofe;i++){
00347 gm[0][i]+=b11[i];
00348 }
00349
00350
00351 vxv (b11,b11,am);
00352 vxm (r,am,av);
00353 for (i=0;i<ndofe;i++){
00354 gm[0][i]+=av[i];
00355 }
00356
00357
00358 vxv (b21,b21,am);
00359 vxm (r,am,av);
00360 for (i=0;i<ndofe;i++){
00361 gm[0][i]+=av[i];
00362 }
00363
00364
00365
00366
00367
00368
00369
00370 for (i=0;i<ndofe;i++){
00371 gm[1][i]+=b22[i];
00372 }
00373
00374
00375 vxv (b22,b22,am);
00376 vxm (r,am,av);
00377 for (i=0;i<ndofe;i++){
00378 gm[1][i]+=av[i];
00379 }
00380
00381
00382 vxv (b12,b12,am);
00383 vxm (r,am,av);
00384 for (i=0;i<ndofe;i++){
00385 gm[1][i]+=av[i];
00386 }
00387
00388
00389
00390
00391
00392
00393
00394 for (i=0;i<ndofe;i++){
00395 gm[2][i]+=b12[i]+b21[i];
00396 }
00397
00398
00399 vxv (b11,b12,am);
00400 vxm (r,am,av);
00401 for (i=0;i<ndofe;i++){
00402 gm[2][i]+=av[i];
00403 }
00404
00405
00406 vxv (b12,b11,am);
00407 vxm (r,am,av);
00408 for (i=0;i<ndofe;i++){
00409 gm[2][i]+=av[i];
00410 }
00411
00412
00413 vxv (b22,b21,am);
00414 vxm (r,am,av);
00415 for (i=0;i<ndofe;i++){
00416 gm[2][i]+=av[i];
00417 }
00418
00419
00420 vxv (b21,b22,am);
00421 vxm (r,am,av);
00422 for (i=0;i<ndofe;i++){
00423 gm[2][i]+=av[i];
00424 }
00425
00426
00427
00428 }
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441 void planeelemlq::gnl_grmatrix (matrix &grm,vector &x,vector &y,double xi,double eta,double &jac)
00442 {
00443 long i;
00444 vector b11(ndofe),b12(ndofe),b21(ndofe),b22(ndofe);
00445
00446 bvectors (x,y,xi,eta,jac,b11,b12,b21,b22);
00447
00448 for (i=0;i<ndofe;i++){
00449 grm[0][i]=b11[i];
00450 grm[1][i]=b12[i];
00451 grm[2][i]=b21[i];
00452 grm[3][i]=b22[i];
00453 }
00454 }
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470 void planeelemlq::dmatblock (long ri,long ci,matrix &d, matrix &dd)
00471 {
00472 fillm (0.0,dd);
00473
00474 if (ri==0 && ci==0){
00475 dd[0][0]=d[0][0]; dd[0][1]=d[0][1];
00476 dd[1][0]=d[1][0]; dd[1][1]=d[1][1];
00477 }
00478 if (ri==0 && ci==1){
00479 dd[0][0]=d[0][2];
00480 dd[1][0]=d[1][2];
00481 }
00482 if (ri==1 && ci==0){
00483 dd[0][0]=d[2][0]; dd[0][1]=d[2][1];
00484 }
00485 if (ri==1 && ci==1){
00486 dd[0][0]=d[2][2];
00487 }
00488 }
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500 void planeelemlq::transf_matrix (ivector &nod,matrix &tmat)
00501 {
00502 long i,n,m;
00503
00504 fillm (0.0,tmat);
00505
00506 n=nod.n;
00507 m=tmat.m;
00508 for (i=0;i<m;i++){
00509 tmat[i][i]=1.0;
00510 }
00511
00512 for (i=0;i<n;i++){
00513 if (Mt->nodes[nod[i]].transf>0){
00514 tmat[i*2][i*2] = Mt->nodes[nod[i]].e1[0]; tmat[i*2][i*2+1] = Mt->nodes[nod[i]].e2[0];
00515 tmat[i*2+1][i*2] = Mt->nodes[nod[i]].e1[1]; tmat[i*2+1][i*2+1] = Mt->nodes[nod[i]].e2[1];
00516 }
00517 }
00518 }
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537 void planeelemlq::gl_stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
00538 {
00539 long i,j,ii,jj,ipp;
00540 double xi,eta,jac,thick;
00541 ivector nodes(nne);
00542 vector w,gp,t(nne);
00543 matrix gmr,gmc,dd,d(tncomp,tncomp);
00544
00545 Mt->give_elemnodes (eid,nodes);
00546 Mc->give_thickness (eid,nodes,t);
00547
00548 fillm (0.0,sm);
00549
00550 for (ii=0;ii<nb;ii++){
00551 allocm (ncomp[ii],ndofe,gmr);
00552 for (jj=0;jj<nb;jj++){
00553 if (intordsm[ii][jj]==0) continue;
00554
00555 allocv (intordsm[ii][jj],w);
00556 allocv (intordsm[ii][jj],gp);
00557
00558 allocm (ncomp[jj],ndofe,gmc);
00559 allocm (ncomp[ii],ncomp[jj],dd);
00560
00561 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00562
00563 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00564
00565 for (i=0;i<intordsm[ii][jj];i++){
00566 xi=gp[i];
00567 for (j=0;j<intordsm[ii][jj];j++){
00568 eta=gp[j];
00569
00570
00571 geom_matrix_block (gmr,ii,x,y,xi,eta,jac);
00572 geom_matrix_block (gmc,jj,x,y,xi,eta,jac);
00573
00574
00575 Mm->matstiff (d,ipp);
00576
00577
00578 dmatblock (ii,jj,d,dd);
00579
00580
00581 thick = approx (xi,eta,t);
00582
00583 jac*=thick*w[i]*w[j];
00584
00585 jac=fabs(jac);
00586
00587
00588 bdbjac (sm,gmr,dd,gmc,jac);
00589
00590 ipp++;
00591 }
00592 }
00593
00594 destrm (dd); destrm (gmc); destrv (gp); destrv (w);
00595 }
00596 destrm (gmr);
00597 }
00598
00599 }
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 void planeelemlq::gnl_stiffness_matrix (long lcid,long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
00615 {
00616 long i,j,ipp;
00617 double xi,eta,jac,jac2,thick;
00618 ivector nodes(nne);
00619 vector w,gp,t(nne),sig(tncomp),r(ndofe);
00620 matrix gm(tncomp,ndofe),grm(4,ndofe),d(tncomp,tncomp),s(4,4);
00621
00622
00623 Mt->give_elemnodes (eid,nodes);
00624
00625 Mc->give_thickness (eid,nodes,t);
00626
00627 eldispl (lcid,eid,r.a);
00628
00629
00630 fillm (0.0,sm);
00631
00632
00633 allocv (intordsm[0][0],w);
00634
00635 allocv (intordsm[0][0],gp);
00636
00637
00638 gauss_points (gp.a,w.a,intordsm[0][0]);
00639
00640
00641 ipp=Mt->elements[eid].ipp[ri][ci];
00642
00643
00644 for (i=0;i<intordsm[0][0];i++){
00645 xi=gp[i];
00646 for (j=0;j<intordsm[0][0];j++){
00647 eta=gp[j];
00648
00649
00650
00651
00652
00653
00654 gngeom_matrix (gm,r,x,y,xi,eta,jac);
00655
00656
00657 Mm->matstiff (d,ipp);
00658
00659
00660 thick = approx (xi,eta,t);
00661
00662 jac*=thick*w[i]*w[j];
00663
00664
00665 bdbjac (sm,gm,d,gm,jac);
00666
00667
00668
00669
00670
00671
00672
00673 gnl_grmatrix (grm,x,y,xi,eta,jac2);
00674
00675
00676 Mm->givestress (lcid,ipp,sig);
00677
00678 s[0][0]=sig[0]; s[0][1]=sig[2];
00679 s[1][0]=sig[2]; s[1][1]=sig[1];
00680
00681 s[2][2]=sig[0]; s[2][3]=sig[2];
00682 s[3][2]=sig[2]; s[3][3]=sig[1];
00683
00684
00685 bdbjac (sm,grm,s,grm,jac);
00686
00687 ipp++;
00688 }
00689 }
00690
00691 destrv (gp); destrv (w);
00692
00693 }
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704 void planeelemlq::res_stiffness_matrix (long lcid,long eid,matrix &sm)
00705 {
00706 long transf;
00707 ivector nodes(nne);
00708 vector x(nne),y(nne);
00709
00710 Mt->give_node_coord2d (x,y,eid);
00711
00712
00713 gl_stiffness_matrix (eid,0,0,sm,x,y);
00714
00715
00716
00717
00718 Mt->give_elemnodes (eid,nodes);
00719 transf = Mt->locsystems (nodes);
00720 if (transf>0){
00721 matrix tmat (ndofe,ndofe);
00722 transf_matrix (nodes,tmat);
00723 glmatrixtransf (sm,tmat);
00724 }
00725
00726 }
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741 void planeelemlq::mass_matrix (long eid,matrix &mm,vector &x,vector &y)
00742 {
00743 long i,j;
00744 double jac,xi,eta,thick,rho;
00745 ivector nodes(nne);
00746 vector w(intordmm),gp(intordmm),t(nne),dens(nne);
00747 matrix n(napfun,ndofe);
00748
00749 Mt->give_elemnodes (eid,nodes);
00750 Mc->give_thickness (eid,nodes,t);
00751 Mc->give_density (eid,nodes,dens);
00752 gauss_points (gp.a,w.a,intordmm);
00753
00754 fillm (0.0,mm);
00755
00756 for (i=0;i<intordmm;i++){
00757 xi=gp[i];
00758 for (j=0;j<intordmm;j++){
00759 eta=gp[j];
00760 jac_2d (jac,x,y,xi,eta);
00761 bf_matrix (n,xi,eta);
00762
00763 thick = approx (xi,eta,t);
00764 rho = approx (xi,eta,dens);
00765 jac*=w[i]*w[j]*thick*rho;
00766
00767 nnj (mm.a,n.a,jac,n.m,n.n);
00768 }
00769 }
00770
00771 }
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782 void planeelemlq::res_mass_matrix (long eid,matrix &mm)
00783 {
00784 long transf;
00785 ivector nodes(nne);
00786 vector x(nne),y(nne);
00787
00788 Mt->give_node_coord2d (x,y,eid);
00789
00790 mass_matrix (eid,mm,x,y);
00791
00792
00793
00794 Mt->give_elemnodes (eid,nodes);
00795 transf = Mt->locsystems (nodes);
00796 if (transf>0){
00797 matrix tmat (ndofe,ndofe);
00798 transf_matrix (nodes,tmat);
00799 glmatrixtransf (mm,tmat);
00800 }
00801 }
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816 void planeelemlq::res_ip_strains (long lcid,long eid)
00817 {
00818 vector aux,x(nne),y(nne),r(ndofe);
00819 ivector cn(ndofe),nodes(nne);
00820 matrix tmat;
00821
00822 Mt->give_elemnodes (eid,nodes);
00823 eldispl (lcid,eid,r.a);
00824
00825
00826 long transf = Mt->locsystems (nodes);
00827 if (transf>0){
00828 allocv (ndofe,aux);
00829 allocm (ndofe,ndofe,tmat);
00830 transf_matrix (nodes,tmat);
00831
00832 lgvectortransf (aux,r,tmat);
00833 copyv (aux,r);
00834 destrv (aux);
00835 destrm (tmat);
00836 }
00837
00838 Mt->give_node_coord2d (x,y,eid);
00839 gl_ip_strains (lcid,eid,0,0,x,y,r);
00840
00841
00842 }
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856 void planeelemlq::gl_ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00857 {
00858 long i,j,ii,ipp;
00859 double xi,eta,jac;
00860 vector gp,w,eps,epscum;
00861 matrix gm;
00862
00863
00864 for (ii=0;ii<nb;ii++){
00865
00866 allocv (intordsm[ii][ii],gp);
00867 allocv (intordsm[ii][ii],w);
00868 allocv (ncomp[ii],eps);
00869 allocm (ncomp[ii],ndofe,gm);
00870
00871 gauss_points (gp.a,w.a,intordsm[ii][ii]);
00872
00873 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00874 for (i=0;i<intordsm[ii][ii];i++){
00875 xi=gp[i];
00876 for (j=0;j<intordsm[ii][ii];j++){
00877 eta=gp[j];
00878
00879 geom_matrix_block (gm,ii,x,y,xi,eta,jac);
00880 mxv (gm,r,eps);
00881
00882 Mm->storestrain (lcid,ipp,cncomp[ii],eps);
00883 ipp++;
00884 }
00885 }
00886
00887 destrm (gm); destrv (eps); destrv (w); destrv (gp);
00888 }
00889
00890
00891
00892
00893
00894 ipp=Mt->elements[eid].ipp[ri][ci];
00895 allocv (ncomp[0],epscum);
00896 allocv (ncomp[0],eps);
00897 for (i=0;i<intordsm[0][0];i++){
00898 for (j=0;j<intordsm[0][0];j++){
00899 Mm->givestrain (lcid,ipp,cncomp[0],eps);
00900 addv (eps,epscum,epscum);
00901 ipp++;
00902 }
00903 }
00904
00905 cmulv (1.0/4.0,epscum,epscum);
00906 ipp=Mt->elements[eid].ipp[ri+1][ci+1];
00907 Mm->storestrain (lcid,ipp,epscum);
00908 destrv (eps);
00909 destrv (epscum);
00910
00911
00912 allocv (ncomp[1],eps);
00913 Mm->givestrain (lcid,ipp,cncomp[1],eps);
00914 ipp=Mt->elements[eid].ipp[ri][ci];
00915 Mm->storestrain (lcid,ipp,cncomp[1],eps); ipp++;
00916 Mm->storestrain (lcid,ipp,cncomp[1],eps); ipp++;
00917 Mm->storestrain (lcid,ipp,cncomp[1],eps); ipp++;
00918 Mm->storestrain (lcid,ipp,cncomp[1],eps); ipp++;
00919 destrv (eps);
00920 }
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934 void planeelemlq::gnl_ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00935 {
00936 long i,j,ipp;
00937 double xi,eta,jac,b11r,b12r,b21r,b22r;
00938 vector gp,w,eps,b11(ndofe),b12(ndofe),b21(ndofe),b22(ndofe);
00939
00940 allocv (intordsm[0][0],gp);
00941 allocv (intordsm[0][0],w);
00942 allocv (tncomp,eps);
00943
00944 gauss_points (gp.a,w.a,intordsm[0][0]);
00945
00946 ipp=Mt->elements[eid].ipp[ri][ci];
00947 for (i=0;i<intordsm[0][0];i++){
00948 xi=gp[i];
00949 for (j=0;j<intordsm[0][0];j++){
00950 eta=gp[j];
00951
00952 bvectors (x,y,xi,eta,jac,b11,b12,b21,b22);
00953
00954 scprd (b11,r,b11r);
00955 scprd (b12,r,b12r);
00956 scprd (b21,r,b21r);
00957 scprd (b22,r,b22r);
00958
00959 eps[0]=b11r+0.5*b11r*b11r+0.5*b21r*b21r;
00960 eps[1]=b22r+0.5*b12r*b12r+0.5*b22r*b22r;
00961 eps[2]=b12r+b21r+b11r*b12r+b21r*b22r;
00962
00963
00964 Mm->storestrain (lcid,ipp,eps);
00965 ipp++;
00966 }
00967 }
00968
00969 destrv (eps); destrv (w); destrv (gp);
00970
00971 }
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983 void planeelemlq::nod_strains_ip (long lcid,long eid,long ri,long ci)
00984 {
00985 long i,j,ipp;
00986 ivector ipnum(nne),nod(nne);
00987 vector eps(gncomp);
00988
00989
00990
00991 ipp=Mt->elements[eid].ipp[ri][ci];
00992 nodip_planelq (ipp,intordsm[0][0],ipnum);
00993
00994
00995 Mt->give_elemnodes (eid,nod);
00996
00997 for (i=0;i<nne;i++){
00998
00999 Mm->givestrain (lcid,ipnum[i],eps);
01000
01001
01002 j=nod[i];
01003 Mt->nodes[j].storestrain (lcid,0,eps);
01004 }
01005
01006 }
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016 void planeelemlq::nod_strains_comp (long lcid,long eid)
01017 {
01018 long i;
01019 double jac;
01020 ivector enod(nne);
01021 vector x(nne),y(nne),nxi(nne),neta(nne),r(ndofe),eps(tncomp),aux;
01022 matrix tmat,gm(tncomp,ndofe);
01023
01024
01025 Mt->give_node_coord2d (x, y, eid);
01026
01027 Mt->give_elemnodes (eid, enod);
01028
01029 eldispl (lcid,eid,r.a);
01030
01031
01032 long transf = Mt->locsystems (enod);
01033 if (transf>0){
01034 allocv (ndofe,aux);
01035 allocm (ndofe,ndofe,tmat);
01036 transf_matrix (enod,tmat);
01037
01038 lgvectortransf (aux,r,tmat);
01039 copyv (aux,r);
01040 destrv (aux);
01041 destrm (tmat);
01042 }
01043
01044
01045
01046 nodcoord_planelq (nxi,neta);
01047
01048
01049 for (i=0;i<nne;i++){
01050
01051 geom_matrix (gm,x,y,nxi[i],neta[i],jac);
01052
01053 mxv (gm,r,eps);
01054
01055
01056 Mt->nodes[enod[i]].storestrain (lcid,0,eps);
01057 }
01058
01059 }
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072 void planeelemlq::strains (long lcid,long eid,long ri,long ci)
01073 {
01074 long i,naep,ncp,sid;
01075 vector coord,eps;
01076
01077 switch (Mm->stra.tape[eid]){
01078 case nowhere:{
01079 break;
01080 }
01081 case intpts:{
01082
01083 res_ip_strains (lcid,eid);
01084 break;
01085 }
01086 case enodes:{
01087
01088 nod_strains_ip (lcid,eid,ri,ci);
01089 break;
01090 }
01091 case cenodes:{
01092
01093 nod_strains_comp (lcid,eid);
01094 break;
01095 }
01096 case userdefined:{
01097
01098 naep = Mm->stra.give_naep (eid);
01099 ncp = Mm->stra.give_ncomp (eid);
01100 sid = Mm->stra.give_sid (eid);
01101 allocv (ncp,eps);
01102 allocv (2,coord);
01103 for (i=0;i<naep;i++){
01104 Mm->stra.give_aepcoord (sid,i,coord);
01105
01106 if (Mp->strainaver==0)
01107
01108 if (Mp->strainaver==1)
01109
01110
01111 Mm->stra.storevalues(lcid,eid,i,eps);
01112 }
01113 destrv (eps);
01114 destrv (coord);
01115 break;
01116 }
01117 default:{
01118 print_err("unknown strain point is required", __FILE__, __LINE__, __func__);
01119 }
01120 }
01121
01122 }
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133 void planeelemlq::res_ip_stresses (long lcid,long eid)
01134 {
01135 long ri,ci;
01136 ri=0;
01137 ci=0;
01138 compute_nlstress (lcid,eid,ri,ci);
01139 }
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151 void planeelemlq::nod_stresses_ip (long lcid,long eid,long ri,long ci)
01152 {
01153 long i,j,ipp;
01154 ivector ipnum(nne),nod(nne);
01155 vector sig(gncomp);
01156
01157
01158
01159 ipp=Mt->elements[eid].ipp[ri][ci];
01160 nodip_planelq (ipp,intordsm[0][0],ipnum);
01161
01162
01163 Mt->give_elemnodes (eid,nod);
01164
01165 for (i=0;i<nne;i++){
01166
01167 Mm->givestress (lcid,ipnum[i],sig);
01168
01169
01170 j=nod[i];
01171 Mt->nodes[j].storestress (lcid,0,sig);
01172 }
01173
01174 }
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185 void planeelemlq::nod_stresses_comp (long lcid,long eid,long ri,long ci,double **stra,double **stre)
01186 {
01187 long i,j,ipp;
01188 vector eps(tncomp),sig(tncomp);
01189 matrix d(tncomp,tncomp);
01190
01191
01192 ipp=Mt->elements[eid].ipp[ri][ci];
01193
01194 Mm->matstiff (d,ipp);
01195
01196
01197 for (i=0;i<nne;i++){
01198 for (j=0;j<eps.n;j++){
01199 eps[j]=stra[i][j];
01200 }
01201 mxv (d,eps,sig);
01202 for (j=0;j<eps.n;j++){
01203 stre[i][j]=sig[j];
01204 }
01205 }
01206 }
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216 void planeelemlq::stresses (long lcid,long eid,long ri,long ci)
01217 {
01218 long i,naep,ncp,sid;
01219 vector coord,sig;
01220 double **stra,**stre;
01221
01222 if (Mp->stressaver==0){
01223 stra = new double* [nne];
01224 stre = new double* [nne];
01225 for (i=0;i<nne;i++){
01226 stra[i] = new double [tncomp];
01227 stre[i] = new double [tncomp];
01228 }
01229
01230
01231 }
01232
01233 res_ip_stresses (lcid,eid);
01234
01235
01236 switch (Mm->stre.tape[eid]){
01237 case nowhere:{
01238 break;
01239 }
01240 case intpts:{
01241 res_ip_stresses (lcid,eid);
01242 break;
01243 }
01244 case enodes:{
01245 nod_stresses_ip (lcid,eid,ri,ci);
01246 break;
01247 }
01248 case userdefined:{
01249
01250 naep = Mm->stre.give_naep (eid);
01251 ncp = Mm->stre.give_ncomp (eid);
01252 sid = Mm->stre.give_sid (eid);
01253 allocv (ncp,sig);
01254 allocv (2,coord);
01255 for (i=0;i<naep;i++){
01256 Mm->stre.give_aepcoord (sid,i,coord);
01257
01258 if (Mp->stressaver==0)
01259
01260 if (Mp->stressaver==1)
01261
01262
01263 Mm->stre.storevalues(lcid,eid,i,sig);
01264 }
01265 destrv (sig);
01266 destrv (coord);
01267 break;
01268 }
01269 default:{
01270 print_err("unknown stress point is required", __FILE__, __LINE__, __func__);
01271 }
01272 }
01273
01274 }
01275
01276
01277
01278
01279
01280
01281
01282
01283 void planeelemlq::nod_eqother_ip (long eid,long ri,long ci)
01284 {
01285 long i,j,ipp,ncompo;
01286 ivector ipnum(nne),nod(nne);
01287 vector eqother;
01288
01289
01290
01291 ipp=Mt->elements[eid].ipp[ri][ci];
01292 nodip_planelq (ipp,intordsm[0][0],ipnum);
01293
01294
01295
01296 Mt->give_elemnodes (eid,nod);
01297
01298 for (i=0;i<nne;i++){
01299 ncompo = Mm->givencompeqother (ipnum[i],0);
01300 reallocv (ncompo,eqother);
01301 Mm->giveeqother (ipnum[i],0,ncompo,eqother.a);
01302
01303
01304 j=nod[i];
01305 Mt->nodes[j].storeother (0,ncompo,eqother);
01306 }
01307 }
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325 void planeelemlq::load_matrix (long eid,matrix &lm,vector &x,vector &y)
01326 {
01327 long i,j;
01328 double jac,xi,eta,w1,w2,thick;
01329 ivector nodes(nne);
01330 vector w(intordmm),gp(intordmm),t(nne);
01331 matrix n(napfun,ndofe);
01332
01333 Mt->give_elemnodes (eid,nodes);
01334 Mc->give_thickness (eid,nodes,t);
01335
01336 gauss_points (gp.a,w.a,intordmm);
01337
01338 fillm (0.0,lm);
01339
01340 for (i=0;i<intordmm;i++){
01341 xi=gp[i]; w1=w[i];
01342 for (j=0;j<intordmm;j++){
01343 eta=gp[j]; w2=w[j];
01344 jac_2d (jac,x,y,xi,eta);
01345 bf_matrix (n,xi,eta);
01346
01347 thick = approx (xi,eta,t);
01348 jac*=w1*w2*thick;
01349
01350 nnj (lm.a,n.a,jac,n.m,n.n);
01351 }
01352 }
01353
01354 }
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367 void planeelemlq::res_load_matrix (long eid,matrix &lm)
01368 {
01369 long transf;
01370 ivector nodes(nne);
01371 vector x(nne),y(nne);
01372
01373 Mt->give_node_coord2d (x,y,eid);
01374
01375 load_matrix (eid,lm,x,y);
01376
01377
01378
01379 Mt->give_elemnodes (eid,nodes);
01380 transf = Mt->locsystems (nodes);
01381 if (transf>0){
01382 matrix tmat (ndofe,ndofe);
01383 transf_matrix (nodes,tmat);
01384 glmatrixtransf (lm,tmat);
01385 }
01386
01387 }
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412 void planeelemlq::gl_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01413 {
01414 integratedquant iq;
01415 iq=locstress;
01416
01417
01418 compute_nlstress (lcid,eid,ri,ci);
01419
01420
01421 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01422 }
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437 void planeelemlq::gnl_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01438 {
01439 long i,j,k,ipp;
01440 double xi,eta,jac,thick;
01441 ivector nodes(nne);
01442 vector w,gp,t(nne),sig(tncomp),contr(ndofe),r(ndofe);
01443 matrix gm(tncomp,ndofe);
01444
01445
01446 Mt->give_elemnodes (eid,nodes);
01447
01448 Mc->give_thickness (eid,nodes,t);
01449
01450 eldispl (lcid,eid,r.a);
01451
01452
01453 fillv (0.0,ifor);
01454
01455
01456 allocv (intordsm[0][0],gp);
01457
01458 allocv (intordsm[0][0],w);
01459
01460
01461 gauss_points (gp.a,w.a,intordsm[0][0]);
01462
01463
01464 ipp=Mt->elements[eid].ipp[ri][ci];
01465
01466
01467 for (i=0;i<intordsm[0][0];i++){
01468 xi=gp[i];
01469 for (j=0;j<intordsm[0][0];j++){
01470 eta=gp[j];
01471
01472
01473 thick = approx (xi,eta,t);
01474
01475
01476 if (Mp->strcomp==1)
01477 Mm->computenlstresses (ipp);
01478
01479 Mm->givestress (lcid,ipp,sig);
01480
01481
01482 gngeom_matrix (gm,r,x,y,xi,eta,jac);
01483
01484 mtxv (gm,sig,contr);
01485
01486 cmulv (jac*w[i]*w[j]*thick,contr);
01487
01488 for (k=0;k<contr.n;k++){
01489 ifor[k]+=contr[k];
01490 }
01491
01492 ipp++;
01493 }
01494 }
01495
01496 destrv (w); destrv (gp);
01497 }
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515 void planeelemlq::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01516 {
01517 integratedquant iq;
01518 iq=nonlocstress;
01519
01520
01521 compute_nonloc_nlstress (lcid,eid,ri,ci);
01522
01523
01524 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01525 }
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540 void planeelemlq::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01541 {
01542 integratedquant iq;
01543 iq=stressincr;
01544
01545
01546 compute_nlstressincr (lcid,eid,ri,ci);
01547
01548
01549 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01550 }
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564 void planeelemlq::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y)
01565 {
01566 integratedquant iq;
01567 iq=eigstress;
01568
01569
01570 if (Mp->eigstrcomp)
01571 compute_eigstress (lcid,eid,ri,ci);
01572
01573
01574 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y);
01575 }
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588 void planeelemlq::res_internal_forces (long lcid,long eid,vector &ifor)
01589 {
01590 long transf;
01591 ivector nodes (nne);
01592 vector v(ndofe),x(nne),y(nne);
01593
01594 Mt->give_node_coord2d (x,y,eid);
01595
01596 gl_internal_forces (lcid,eid,0,0,ifor,x,y);
01597
01598
01599
01600 Mt->give_elemnodes (eid,nodes);
01601 transf = Mt->locsystems (nodes);
01602 if (transf>0){
01603 matrix tmat (ndofe,ndofe);
01604 transf_matrix (nodes,tmat);
01605
01606 glvectortransf (ifor,v,tmat);
01607 copyv (v,ifor);
01608 }
01609 }
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622 void planeelemlq::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
01623 {
01624 long transf;
01625 ivector nodes (nne);
01626 vector v(ndofe),x(nne),y(nne);
01627
01628 Mt->give_node_coord2d (x,y,eid);
01629
01630 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y);
01631
01632
01633
01634 Mt->give_elemnodes (eid,nodes);
01635 transf = Mt->locsystems (nodes);
01636 if (transf>0){
01637 matrix tmat (ndofe,ndofe);
01638 transf_matrix (nodes,tmat);
01639
01640 glvectortransf (ifor,v,tmat);
01641 copyv (v,ifor);
01642 }
01643 }
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656 void planeelemlq::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01657 {
01658 long transf;
01659 ivector nodes (nne);
01660 vector v(ndofe),x(nne),y(nne);
01661
01662 Mt->give_node_coord2d (x,y,eid);
01663
01664 incr_internal_forces (lcid,eid,0,0,ifor,x,y);
01665
01666
01667
01668 Mt->give_elemnodes (eid,nodes);
01669 transf = Mt->locsystems (nodes);
01670 if (transf>0){
01671 matrix tmat (ndofe,ndofe);
01672 transf_matrix (nodes,tmat);
01673
01674 glvectortransf (ifor,v,tmat);
01675 copyv (v,ifor);
01676 }
01677 }
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690 void planeelemlq::res_eigstrain_forces (long lcid,long eid,vector &nfor)
01691 {
01692 long transf;
01693 ivector nodes (nne);
01694 vector v(ndofe),x(nne),y(nne);
01695
01696 Mt->give_node_coord2d (x,y,eid);
01697
01698 eigstrain_forces (lcid,eid,0,0,nfor,x,y);
01699
01700
01701
01702 Mt->give_elemnodes (eid,nodes);
01703 transf = Mt->locsystems (nodes);
01704 if (transf>0){
01705 matrix tmat (ndofe,ndofe);
01706 transf_matrix (nodes,tmat);
01707
01708 glvectortransf (nfor,v,tmat);
01709 copyv (v,nfor);
01710 }
01711 }
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725 void planeelemlq::local_values(long lcid,long eid,long ri,long ci)
01726 {
01727 long i,j,ii,ipp;
01728
01729 for (ii=0;ii<nb;ii++){
01730 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01731
01732 for (i=0;i<intordsm[ii][ii];i++){
01733 for (j=0;j<intordsm[ii][ii];j++){
01734
01735 if (Mp->strcomp==1)
01736 Mm->computenlstresses (ipp);
01737
01738 ipp++;
01739 }
01740 }
01741 }
01742 }
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755 void planeelemlq::compute_nlstress (long lcid,long eid,long ri,long ci)
01756 {
01757 long i,j,ii,ipp;
01758
01759 for (ii=0;ii<nb;ii++){
01760 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01761
01762 for (i=0;i<intordsm[ii][ii];i++){
01763 for (j=0;j<intordsm[ii][ii];j++){
01764
01765 if (Mp->strcomp==1)
01766 Mm->computenlstresses (ipp);
01767
01768 ipp++;
01769 }
01770 }
01771 }
01772 }
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785 void planeelemlq::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01786 {
01787 long i,j,ii,ipp;
01788
01789 for (ii=0;ii<nb;ii++){
01790 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01791
01792 for (i=0;i<intordsm[ii][ii];i++){
01793 for (j=0;j<intordsm[ii][ii];j++){
01794
01795 if (Mp->strcomp==1)
01796 Mm->computenlstressesincr (ipp);
01797 ipp++;
01798 }
01799 }
01800 }
01801 }
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814 void planeelemlq::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01815 {
01816 long i,j,ii,ipp;
01817
01818 for (ii=0;ii<nb;ii++){
01819 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01820
01821 for (i=0;i<intordsm[ii][ii];i++){
01822 for (j=0;j<intordsm[ii][ii];j++){
01823
01824
01825 if (Mp->strcomp==1)
01826 Mm->compnonloc_nlstresses (ipp);
01827
01828 ipp++;
01829 }
01830 }
01831 }
01832 }
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843 void planeelemlq::compute_eigstress (long lcid,long eid,long ri,long ci)
01844 {
01845 long i,j,ii,ipp;
01846 vector eigstr(tncomp),sig(tncomp);
01847 matrix d(tncomp,tncomp);
01848
01849 for (ii=0;ii<nb;ii++){
01850 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01851
01852 for (i=0;i<intordsm[ii][ii];i++){
01853 for (j=0;j<intordsm[ii][ii];j++){
01854
01855 Mm->giveeigstrain (ipp,eigstr);
01856
01857
01858 Mm->matstiff (d,ipp);
01859
01860 mxv (d,eigstr,sig);
01861
01862 Mm->storeeigstress (ipp,sig);
01863
01864 ipp++;
01865 }
01866 }
01867 }
01868 }
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885 void planeelemlq::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y)
01886 {
01887 long i,j,ii,ipp;
01888 double xi,eta,jac,thick;
01889 ivector nodes(nne);
01890 vector w,gp,t(nne),ipv,contr(ndofe);
01891 matrix gm;
01892
01893 Mt->give_elemnodes (eid,nodes);
01894 Mc->give_thickness (eid,nodes,t);
01895
01896 fillv (0.0,nv);
01897
01898 for (ii=0;ii<nb;ii++){
01899 allocv (intordsm[ii][ii],gp);
01900 allocv (intordsm[ii][ii],w);
01901 allocm (ncomp[ii],ndofe,gm);
01902 allocv (ncomp[ii],ipv);
01903
01904 gauss_points (gp.a,w.a,intordsm[ii][ii]);
01905 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01906
01907 for (i=0;i<intordsm[ii][ii];i++){
01908 xi=gp[i];
01909 for (j=0;j<intordsm[ii][ii];j++){
01910 eta=gp[j];
01911 thick = approx (xi,eta,t);
01912
01913
01914 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
01915
01916
01917 geom_matrix_block (gm,ii,x,y,xi,eta,jac);
01918
01919
01920 mtxv (gm,ipv,contr);
01921
01922 cmulv (jac*w[i]*w[j]*thick,contr);
01923
01924
01925 addv(contr,nv,nv);
01926
01927 ipp++;
01928 }
01929 }
01930 destrm (gm); destrv (ipv); destrv (w); destrv (gp);
01931 }
01932 }
01933
01934
01935 void planeelemlq::nodeforces (long eid,long *le,double *nv,vector &nf)
01936 {
01937 long i;
01938 double ww,jac,xi,eta;
01939 vector x(nne),y(nne),gp(intordb),w(intordb),av(ndofe),v(ndofe);
01940 matrix n(napfun,ndofe),am(ndofe,ndofe);
01941
01942 Mt->give_node_coord2d (x,y,eid);
01943 gauss_points (gp.a,w.a,intordb);
01944
01945 if (le[0]==1){
01946 fillm (0.0,am);
01947 eta = 1.0;
01948 for (i=0;i<intordb;i++){
01949 xi = gp[i];
01950 ww = w[i];
01951
01952 bf_matrix (n,xi,eta);
01953
01954 jac1d_2d (jac,x,y,xi,0);
01955 jac *= ww;
01956
01957 nnj (am.a,n.a,jac,n.m,n.n);
01958 }
01959 fillv (0.0,av);
01960 av[0]=nv[0]; av[1]=nv[1]; av[2]=nv[2]; av[3]=nv[3];
01961 mxv (am,av,v); addv (nf,v,nf);
01962 }
01963 if (le[1]==1){
01964 fillm (0.0,am);
01965 xi = -1.0;
01966 for (i=0;i<intordb;i++){
01967 eta = gp[i];
01968 ww = w[i];
01969
01970 bf_matrix (n,xi,eta);
01971
01972 jac1d_2d (jac,x,y,eta,1);
01973 jac *= ww;
01974
01975 nnj (am.a,n.a,jac,n.m,n.n);
01976 }
01977 fillv (0.0,av);
01978 av[2]=nv[4]; av[3]=nv[5]; av[4]=nv[6]; av[5]=nv[7];
01979 mxv (am,av,v); addv (nf,v,nf);
01980 }
01981 if (le[2]==1){
01982 fillm (0.0,am);
01983 eta = -1.0;
01984 for (i=0;i<intordb;i++){
01985 xi = gp[i];
01986 ww = w[i];
01987
01988 bf_matrix (n,xi,eta);
01989
01990 jac1d_2d (jac,x,y,xi,2);
01991 jac *= ww;
01992
01993 nnj (am.a,n.a,jac,n.m,n.n);
01994 }
01995 fillv (0.0,av);
01996 av[4]=nv[8]; av[5]=nv[9]; av[6]=nv[10]; av[7]=nv[11];
01997 mxv (am,av,v); addv (nf,v,nf);
01998 }
01999 if (le[3]==1){
02000 fillm (0.0,am);
02001 xi = 1.0;
02002 for (i=0;i<intordb;i++){
02003 eta = gp[i];
02004 ww = w[i];
02005
02006 bf_matrix (n,xi,eta);
02007
02008 jac1d_2d (jac,x,y,eta,3);
02009 jac *= ww;
02010
02011 nnj (am.a,n.a,jac,n.m,n.n);
02012 }
02013 fillv (0.0,av);
02014 av[6]=nv[12]; av[7]=nv[13]; av[0]=nv[14]; av[1]=nv[15];
02015 mxv (am,av,v); addv (nf,v,nf);
02016 }
02017 }
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033 void planeelemlq::edge_integral (long edg,vector &x,vector &y,long intord,vector &gp,vector &w,
02034 vector &t,vector &coef,matrix &km)
02035 {
02036 long i;
02037 double xi,eta,jac,ipval,thick;
02038 matrix n(napfun,ndofe);
02039
02040 if (edg==0){
02041 eta=1.0;
02042 for (i=0;i<intord;i++){
02043 xi=gp[i];
02044
02045 bf_matrix (n,xi,eta);
02046 jac1d_2d (jac,x,y,xi,edg);
02047
02048 ipval=approx (xi,eta,coef);
02049 thick=approx (xi,eta,t);
02050
02051 jac*=w[i]*ipval*thick;
02052 nnj (km.a,n.a,jac,n.m,n.n);
02053 }
02054 }
02055
02056 if (edg==1){
02057 xi=-1.0;
02058 for (i=0;i<intord;i++){
02059 eta=gp[i];
02060
02061 bf_matrix (n,xi,eta);
02062 jac1d_2d (jac,x,y,eta,edg);
02063
02064 ipval=approx (xi,eta,coef);
02065 thick=approx (xi,eta,t);
02066
02067 jac*=w[i]*ipval*thick;
02068 nnj (km.a,n.a,jac,n.m,n.n);
02069 }
02070 }
02071
02072 if (edg==2){
02073 eta=-1.0;
02074 for (i=0;i<intord;i++){
02075 xi=gp[i];
02076
02077 bf_matrix (n,xi,eta);
02078 jac1d_2d (jac,x,y,xi,edg);
02079
02080 ipval=approx (xi,eta,coef);
02081 thick=approx (xi,eta,t);
02082
02083 jac*=w[i]*ipval*thick;
02084 nnj (km.a,n.a,jac,n.m,n.n);
02085 }
02086 }
02087
02088 if (edg==3){
02089 xi=1.0;
02090 for (i=0;i<intord;i++){
02091 eta=gp[i];
02092
02093 bf_matrix (n,xi,eta);
02094 jac1d_2d (jac,x,y,eta,edg);
02095
02096 ipval=approx (xi,eta,coef);
02097 thick=approx (xi,eta,t);
02098
02099 jac*=w[i]*ipval*thick;
02100 nnj (km.a,n.a,jac,n.m,n.n);
02101 }
02102 }
02103
02104 }
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118 void planeelemlq::edgenodeval (long edg,vector &nodval,double *list)
02119 {
02120 long i,j,k;
02121 ivector edgenod(nned);
02122
02123 fillv (0.0,nodval);
02124 linquadrilat_edgnod (edgenod.a,edg);
02125
02126 k=0;
02127 for (i=0;i<nned;i++){
02128 for (j=0;j<napfun;j++){
02129 nodval[k]=list[edg*nned*napfun+k];
02130 k++;
02131 }
02132 }
02133 }
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159 void planeelemlq::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
02160 {
02161 long i, j, k, l, ipp;
02162 long ii, jj, nv = nodval.n;
02163 double xi, eta, ipval;
02164 vector w, gp, anv(nne);
02165 long nstra, nstre, ncompstr, ncompeqother;
02166 long idstra, idstre, idoth, idic;
02167 inictype ict;
02168
02169 nstra = idstra = nstre = idstre = idoth = idic = 0;
02170
02171 ict = ictn[0];
02172 for (i=0; i<nne; i++)
02173 {
02174 if (ictn[i] != ict)
02175 {
02176 print_err("Incompatible types of initial conditions on element %ld\n"
02177 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
02178 abort();
02179 }
02180 }
02181 for (j = 0; j < nv; j++)
02182 {
02183 for(i = 0; i < nne; i++)
02184 anv[i] = nodval[i][j];
02185 for (ii = 0; ii < nb; ii++)
02186 {
02187 for (jj = 0; jj < nb; jj++)
02188 {
02189 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
02190 if (intordsm[ii][jj] == 0)
02191 continue;
02192 allocv (intordsm[ii][jj],gp);
02193 allocv (intordsm[ii][jj],w);
02194 gauss_points (gp.a,w.a,intordsm[ii][jj]);
02195 for (k = 0; k < intordsm[ii][jj]; k++)
02196 {
02197 xi=gp[k];
02198 for (l = 0; l < intordsm[ii][jj]; l++)
02199 {
02200 eta=gp[l];
02201
02202 ipval = approx (xi,eta,anv);
02203 ncompstr = Mm->ip[ipp].ncompstr;
02204 ncompeqother = Mm->ip[ipp].ncompeqother;
02205 if ((ictn[0] & inistrain) && (j < ncompstr))
02206 {
02207 Mm->ip[ipp].strain[idstra] += ipval;
02208 ipp++;
02209 continue;
02210 }
02211 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02212 {
02213 Mm->ip[ipp].stress[idstre] += ipval;
02214 ipp++;
02215 continue;
02216 }
02217 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
02218 {
02219 Mm->ip[ipp].eqother[idoth] += ipval;
02220 ipp++;
02221 continue;
02222 }
02223 if ((ictn[0] & inicond) && (j < nv))
02224 {
02225 if (Mm->ic[ipp] == NULL)
02226 {
02227 Mm->ic[ipp] = new double[nv-j];
02228 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
02229 }
02230 Mm->ic[ipp][idic] += ipval;
02231 ipp++;
02232 continue;
02233 }
02234 ipp++;
02235
02236 }
02237 }
02238 destrv (gp); destrv (w);
02239 }
02240 }
02241 ipp=Mt->elements[eid].ipp[ri][ci];
02242 ncompstr = Mm->ip[ipp].ncompstr;
02243 ncompeqother = Mm->ip[ipp].ncompeqother;
02244 if ((ictn[0] & inistrain) && (j < ncompstr))
02245 {
02246 nstra++;
02247 idstra++;
02248 continue;
02249 }
02250 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02251 {
02252 nstre++;
02253 idstre++;
02254 continue;
02255 }
02256 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
02257 {
02258 idoth++;
02259 continue;
02260 }
02261 if ((ictn[0] & inicond) && (j < nv))
02262 {
02263 idic++;
02264 continue;
02265 }
02266 }
02267 }
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280 void planeelemlq::ipvolume (long eid,long ri,long ci)
02281 {
02282 long i,j,ii,ipp;
02283 double xi,eta,jac,thick;
02284 ivector nodes(nne);
02285 vector x(nne),y(nne),t(nne),w,gp;
02286
02287 Mt->give_node_coord2d (x,y,eid);
02288 Mt->give_elemnodes (eid,nodes);
02289 Mc->give_thickness (eid,nodes,t);
02290
02291 for (ii=0;ii<nb;ii++){
02292 allocv (intordsm[ii][ii],w);
02293 allocv (intordsm[ii][ii],gp);
02294
02295 gauss_points (gp.a,w.a,intordsm[ii][ii]);
02296
02297 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
02298
02299 for (i=0;i<intordsm[ii][ii];i++){
02300 xi=gp[i];
02301 for (j=0;j<intordsm[ii][ii];j++){
02302 eta=gp[j];
02303
02304 thick = approx (xi,eta,t);
02305 jac_2d (jac,x,y,xi,eta);
02306 jac*=w[i]*w[j]*thick;
02307
02308 Mm->storeipvol (ipp,jac);
02309
02310 ipp++;
02311 }
02312 }
02313 destrv (gp); destrv (w);
02314 }
02315 }
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326 void planeelemlq::intpointval (long eid,vector &nodval,vector &ipval)
02327 {
02328 long i,j,ii,k;
02329 double xi,eta;
02330 vector w,gp;
02331
02332 k=0;
02333 for (ii=0;ii<nb;ii++){
02334 allocv (intordsm[ii][ii],w);
02335 allocv (intordsm[ii][ii],gp);
02336
02337 gauss_points (gp.a,w.a,intordsm[ii][ii]);
02338
02339 for (i=0;i<intordsm[ii][ii];i++){
02340 xi=gp[i];
02341 for (j=0;j<intordsm[ii][ii];j++){
02342 eta=gp[j];
02343
02344 ipval[k]=approx (xi,eta,nodval);
02345 k++;
02346 }
02347 }
02348
02349 destrv (w);
02350 destrv (gp);
02351 }
02352 }
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365 void planeelemlq::define_meaning (long eid)
02366 {
02367 ivector cn(ndofe),nod(nne);
02368
02369 Mt->give_elemnodes (eid,nod);
02370 Mt->give_code_numbers (eid,cn.a);
02371
02372
02373 if (cn[0]>0) Mt->nodes[nod[0]].meaning[0]=1;
02374
02375 if (cn[1]>0) Mt->nodes[nod[0]].meaning[1]=2;
02376
02377 if (cn[2]>0) Mt->nodes[nod[1]].meaning[0]=1;
02378
02379 if (cn[3]>0) Mt->nodes[nod[1]].meaning[1]=2;
02380
02381 if (cn[4]>0) Mt->nodes[nod[2]].meaning[0]=1;
02382
02383 if (cn[5]>0) Mt->nodes[nod[2]].meaning[1]=2;
02384
02385 if (cn[6]>0) Mt->nodes[nod[3]].meaning[0]=1;
02386
02387 if (cn[7]>0) Mt->nodes[nod[3]].meaning[1]=2;
02388
02389 }
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411 void planeelemlq::ntdbr_vector (long eid,vector &ntdbr)
02412 {
02413 long intord = 2;
02414 long i,j,k,l,ipp,ippsh,ri,ci,lcid,body;
02415 double thick,xi,eta,jac;
02416 ivector nodes(nne);
02417 vector x(nne),y(nne),gp(intord),w(intord),t(nne),eps(tncomp),sig(tncomp),bf(nne),r;
02418 matrix d(tncomp,tncomp),gm;
02419
02420 ri = ci = lcid = 0;
02421 ipp = Mt->elements[eid].ipp[ri][ci];
02422
02423 body = 1;
02424
02425 Mt->give_elemnodes (eid,nodes);
02426 Mt->give_node_coord2d (x,y,eid);
02427 Mc->give_thickness (eid,nodes,t);
02428
02429 Mm->matstiff (d,ipp);
02430
02431
02432 ippsh = ipp+4;
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448 gauss_points (gp.a,w.a,intord);
02449
02450 fillv (0.0,ntdbr);
02451
02452 for (i=0;i<intord;i++){
02453 xi=gp[i];
02454 for (j=0;j<intord;j++){
02455 eta=gp[j];
02456
02457 eps[0] = Mm->ip[ipp ].strain[0];
02458 eps[1] = Mm->ip[ipp ].strain[1];
02459 eps[2] = Mm->ip[ippsh].strain[2];
02460 ipp++;
02461 jac_2d (jac,x,y,xi,eta);
02462
02463 mxv (d,eps,sig);
02464
02465 bf_lin_4_2d (bf.a,xi,eta);
02466
02467 thick = approx (xi,eta,t);
02468 jac*=thick*w[i]*w[j];
02469
02470 for (k=0;k<tncomp;k++)
02471 for (l=0;l<nne;l++)
02472 ntdbr[k*nne+l] += jac * bf[l] * sig[k];
02473
02474 }
02475 }
02476 }
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488 void planeelemlq::ntn_matrix (long eid,matrix &ntn)
02489 {
02490 long intord = 2;
02491 long i,j,k,l;
02492 double thick,xi,eta,jac;
02493 ivector nodes(nne);
02494 vector x(nne),y(nne),t(nne),gp(intord),w(intord),bf(nne);
02495
02496 Mt->give_elemnodes (eid,nodes);
02497 Mt->give_node_coord2d (x,y,eid);
02498 Mc->give_thickness (eid,nodes,t);
02499
02500 gauss_points (gp.a,w.a,intord);
02501
02502 fillm (0.0,ntn);
02503
02504 for (i=0;i<intord;i++){
02505 xi=gp[i];
02506 for (j=0;j<intord;j++){
02507 eta=gp[j];
02508
02509 thick = approx (xi,eta,t);
02510 jac_2d (jac,x,y,xi,eta);
02511 jac*=thick*w[i]*w[j];
02512
02513 bf_lin_4_2d (bf.a,xi,eta);
02514
02515 for (k=0;k<nne;k++)
02516 for (l=0;l<nne;l++)
02517 ntn[k][l] += jac * bf[k] * bf[l];
02518
02519 }
02520 }
02521 }
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546 double planeelemlq :: compute_error (long eid, double &e2, double &u2, double &sizel, vector *rsigfull)
02547 {
02548 long intord = 2;
02549 long i,j,ipp,ippsh,ri,ci;
02550 double thick,area,xi,eta,jac,contr;
02551 double zero=1.0e-20;
02552 ivector nodes(nne);
02553 vector x(nne),y(nne),t(nne),gp(intord),w(intord),bf(nne),r;
02554 vector sig_star(tncomp),sig_roof(tncomp),sig_err(tncomp),eps(tncomp);
02555 matrix d(tncomp,tncomp),dinv(tncomp,tncomp),gm;
02556
02557 ri = ci = 0;
02558 ipp = Mt->elements[eid].ipp[ri][ci];
02559
02560
02561
02562 Mt->give_elemnodes (eid,nodes);
02563 Mt->give_node_coord2d (x,y,eid);
02564 Mc->give_thickness (eid,nodes,t);
02565
02566 Mm->matstiff (d,ipp);
02567 invm (d,dinv,zero);
02568
02569
02570 ippsh = ipp+4;
02571
02572
02573
02574
02575
02576
02577
02578
02579
02580
02581
02582
02583
02584
02585
02586 gauss_points (gp.a,w.a,intord);
02587
02588 e2 = u2 = 0;
02589
02590 for (i=0;i<intord;i++){
02591 xi=gp[i];
02592 for (j=0;j<intord;j++){
02593 eta=gp[j];
02594
02595
02596 eps[0] = Mm->ip[ipp ].strain[0];
02597 eps[1] = Mm->ip[ipp ].strain[1];
02598 eps[2] = Mm->ip[ippsh].strain[2];
02599 ipp++;
02600 jac_2d (jac,x,y,xi,eta);
02601
02602
02603
02604
02605
02606
02607 mxv (d,eps,sig_roof);
02608
02609 bf_lin_4_2d (bf.a,xi,eta);
02610 give_der_star (bf,rsigfull,nodes,sig_star,Mt->nn);
02611
02612 subv (sig_star,sig_roof,sig_err);
02613
02614 thick = approx (xi,eta,t);
02615 jac*=thick*w[i]*w[j];
02616
02617 vxmxv (sig_err,dinv,contr);
02618 e2 += jac*contr;
02619
02620 vxmxv (eps,d,contr);
02621 u2 += jac*contr;
02622 }
02623 }
02624
02625 area = ( (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]) + (x[2]-x[0])*(y[3]-y[0])-(x[3]-x[0])*(y[2]-y[0]) ) / 2.0;
02626 sizel = sqrt(area);
02627 return area;
02628 }
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640
02641 void planeelemlq :: elchar (long eid, matrix &spsig)
02642 {
02643 long ipp,ri,ci,lcid;
02644 double jac;
02645 vector x(nne),y(nne),r(ndofe),eps(tncomp);
02646 matrix gm(tncomp,ndofe),d(tncomp,tncomp);
02647
02648 allocm(1, tncomp, spsig);
02649
02650 lcid = ri = ci = 0;
02651 ipp = Mt->elements[eid].ipp[ri][ci];
02652
02653 Mm->matstiff (d,ipp);
02654
02655 Mt->give_node_coord2d (x,y,eid);
02656 eldispl (lcid,eid,r.a);
02657
02658 geom_matrix (gm,x,y,0.0,0.0,jac);
02659 mxv (gm,r,eps);
02660 mxv (d.a, eps.a, spsig.a, d.m, d.n);
02661
02662 vector aux(tncomp);
02663 Mm->givestress (0,ipp,aux);
02664 ;
02665 vector aux2(tncomp);
02666 Mm->givestrain (0,ipp,aux2);
02667 ;
02668 }
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681 void planeelemlq::extract (atsel &at,vector &val)
02682 {
02683 long i;
02684
02685 for (i=0;i<at.num;i++){
02686 switch (at.atrib[i]){
02687 case 0:{
02688 break;
02689 }
02690 default:{
02691 fprintf (stderr,"\n\n wrong number of atribute in function extract (file %s, line %d).\n",__FILE__,__LINE__);
02692 }
02693 }
02694 }
02695
02696 }
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710 void planeelemlq::mechq_nodval (long eid,vector &nodval,nontransquant qt)
02711 {
02712 long i,ipid;
02713 ivector ipnum(nne);
02714
02715
02716
02717 ipid=Mt->elements[eid].ipp[0][0];
02718 nodip_planelq (ipid,intordsm[0][0],ipnum);
02719
02720 for (i=0;i<nne;i++){
02721
02722 nodval[i] = Mm->givemechq(qt, ipnum[i]);
02723 }
02724 }
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734
02735
02736
02737
02738
02739
02740
02741
02742
02743 void planeelemlq::mechq_nodval_comp (long eid, vector &nodval, long ncnv, long nq, nontransquant *qt)
02744 {
02745 long i, j, ncompstr, ncompo, ncompnl, ipid;
02746 ivector ipnum(nne), enod(nne);
02747 intpoints ipb;
02748
02749
02750 ipid=Mt->elements[eid].ipp[0][0];
02751
02752
02753 Mt->give_elemnodes(eid, enod);
02754
02755
02756
02757 nodip_planelq (ipid,intordsm[0][0],ipnum);
02758
02759
02760 nod_strains_comp(0, eid);
02761
02762
02763 ncompnl = Mm->give_num_averq(ipid, Mm->givenonlocid(ipid));
02764
02765
02766
02767 ipb.copy(Mm->ip[ipid], Mb->nlc, ncompnl);
02768
02769
02770
02771
02772 for (i=0;i<ncnv;i++)
02773 {
02774
02775 ncompstr = Mm->ip[ipnum[i]].ncompstr;
02776
02777 Mm->storestrain(0, ipid, 0, ncompstr, Mt->nodes[enod[i]].strain);
02778
02779 ncompo = Mm->ip[ipnum[i]].ncompeqother;
02780 if (ncompo){
02781
02782 Mm->storeeqother(ipid, 0, ncompo, Mm->ip[ipnum[i]].eqother);
02783 }
02784
02785 if (ncompnl){
02786
02787 Mm->storenonloc(ipid, 0, ncompnl, Mm->ip[ipnum[i]].nonloc);
02788 }
02789
02790
02791 switch(Mp->matmodel)
02792 {
02793 case local:
02794 Mm->computenlstresses(ipid);
02795 break;
02796 case nonlocal:
02797 Mm->compnonloc_nlstresses (ipid);
02798 break;
02799 default:
02800 print_err("unknown approach in material model is required (Mp->matmodel=%ld)", __FILE__, __LINE__, __func__, Mp->matmodel);
02801 }
02802
02803
02804
02805 Mm->updateipvalmat(ipid, 0, 0);
02806
02807
02808 for (j=0; j<nq; j++)
02809 nodval[j*ncnv+i] = Mm->givemechq(qt[j], ipid);
02810 }
02811
02812
02813 Mm->storestrain(0, ipid, 0, Mb->nlc*ipb.ncompstr, ipb.strain);
02814 Mm->storestress(0, ipid, 0, Mb->nlc*ipb.ncompstr, ipb.stress);
02815 Mm->storeother(ipid, 0, ipb.ncompother, ipb.other);
02816 Mm->storeeqother(ipid, 0, ipb.ncompeqother, ipb.eqother);
02817 Mm->storenonloc(ipid, 0, ncompnl, ipb.nonloc);
02818 }