00001 #include "plelemrotlt.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
00010 planeelemrotlt::planeelemrotlt (void)
00011 {
00012 long i,j;
00013
00014
00015 nne=3;
00016
00017 ndofe=9;
00018
00019 tncomp=3;
00020
00021 gncomp=4;
00022
00023 napfun=3;
00024
00025 intordmm=3;
00026
00027 ned=3;
00028
00029 nned=2;
00030
00031
00032 nb=2;
00033
00034
00035 ncomp = new long [nb];
00036 ncomp[0]=2;
00037 ncomp[1]=1;
00038
00039
00040 cncomp = new long [nb];
00041 cncomp[0]=0;
00042 cncomp[1]=2;
00043
00044
00045
00046 nip = new long* [nb];
00047 intordsm = new long* [nb];
00048 for (i=0;i<nb;i++){
00049 nip[i] = new long [nb];
00050 intordsm[i] = new long [nb];
00051 }
00052
00053 nip[0][0]=3; nip[0][1]=0;
00054 nip[1][0]=0; nip[1][1]=3;
00055
00056
00057 tnip=0;
00058 for (i=0;i<nb;i++){
00059 for (j=0;j<nb;j++){
00060 tnip+=nip[i][j];
00061 }
00062 }
00063
00064 intordsm[0][0]=3; intordsm[0][1]=0;
00065 intordsm[1][0]=0; intordsm[1][1]=3;
00066
00067 }
00068
00069 planeelemrotlt::~planeelemrotlt (void)
00070 {
00071 long i;
00072
00073 for (i=0;i<nb;i++){
00074 delete [] nip[i];
00075 delete [] intordsm[i];
00076 }
00077 delete [] nip;
00078 delete [] intordsm;
00079
00080 delete [] cncomp;
00081 delete [] ncomp;
00082 }
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 double planeelemrotlt::approx (vector &areacoord,vector &nodval)
00095 {
00096 double f;
00097 scprd (areacoord,nodval,f);
00098 return f;
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 double planeelemrotlt::approx_nat (double xi,double eta,vector &nodval)
00110 {
00111 double f;
00112 vector areacoord(3);
00113
00114 areacoord[0]=xi;
00115 areacoord[1]=eta;
00116 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00117
00118 scprd (areacoord,nodval,f);
00119 return f;
00120 }
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 void planeelemrotlt::bf_matrix (matrix &n,vector &l,vector &b,vector &c)
00132 {
00133 vector bf(9);
00134
00135 fillm (0.0,n);
00136
00137 bf_rot_3_2d (bf.a,l.a,b.a,c.a);
00138
00139 n[0][0]=bf[0];
00140 n[0][2]=bf[3];
00141 n[0][3]=bf[1];
00142 n[0][5]=bf[4];
00143 n[0][6]=bf[2];
00144 n[0][8]=bf[5];
00145
00146 n[1][1]=bf[0];
00147 n[1][2]=bf[6];
00148 n[1][4]=bf[1];
00149 n[1][5]=bf[7];
00150 n[1][7]=bf[2];
00151 n[1][8]=bf[8];
00152 }
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 void planeelemrotlt::geom_matrix (matrix &gm,vector &x,vector &y,vector &areacoord)
00164 {
00165 double area,det;
00166 vector b(3),c(3),dx(9),dy(9);
00167
00168
00169 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00170 area=det/2.0;
00171 plsb (b.a,y.a,det);
00172 plsc (c.a,x.a,det);
00173
00174 dx_bf_rot_3_2d (dx.a,areacoord.a,b.a,c.a);
00175 dy_bf_rot_3_2d (dy.a,areacoord.a,b.a,c.a);
00176
00177 fillm (0.0,gm);
00178
00179 gm[0][0]=dx[0];
00180 gm[0][2]=dx[3];
00181 gm[0][3]=dx[1];
00182 gm[0][5]=dx[4];
00183 gm[0][6]=dx[2];
00184 gm[0][8]=dx[5];
00185
00186 gm[1][1]=dy[0];
00187 gm[1][2]=dy[6];
00188 gm[1][4]=dy[1];
00189 gm[1][5]=dy[7];
00190 gm[1][7]=dy[2];
00191 gm[1][8]=dy[8];
00192
00193 gm[2][0]=dy[0];
00194 gm[2][1]=dx[0];
00195 gm[2][2]=dy[3]+dx[6];
00196 gm[2][3]=dy[1];
00197 gm[2][4]=dx[1];
00198 gm[2][5]=dy[4]+dx[7];
00199 gm[2][6]=dy[2];
00200 gm[2][7]=dx[2];
00201 gm[2][8]=dy[5]+dx[8];
00202
00203 }
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215 void planeelemrotlt::geom_matrix_block (matrix &gm,long ri,vector &x,vector &y,vector &areacoord)
00216 {
00217 double area,det;
00218 vector b(3),c(3),dx(9),dy(9);
00219
00220
00221 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00222 area=det/2.0;
00223 plsb (b.a,y.a,det);
00224 plsc (c.a,x.a,det);
00225
00226 dx_bf_rot_3_2d (dx.a,areacoord.a,b.a,c.a);
00227 dy_bf_rot_3_2d (dy.a,areacoord.a,b.a,c.a);
00228
00229 fillm (0.0,gm);
00230
00231 if (ri==0){
00232 gm[0][0]=dx[0];
00233 gm[0][2]=dx[3];
00234 gm[0][3]=dx[1];
00235 gm[0][5]=dx[4];
00236 gm[0][6]=dx[2];
00237 gm[0][8]=dx[5];
00238
00239 gm[1][1]=dy[0];
00240 gm[1][2]=dy[6];
00241 gm[1][4]=dy[1];
00242 gm[1][5]=dy[7];
00243 gm[1][7]=dy[2];
00244 gm[1][8]=dy[8];
00245 }
00246
00247 if (ri==1){
00248 gm[0][0]=dy[0];
00249 gm[0][1]=dx[0];
00250 gm[0][2]=dy[3]+dx[6];
00251 gm[0][3]=dy[1];
00252 gm[0][4]=dx[1];
00253 gm[0][5]=dy[4]+dx[7];
00254 gm[0][6]=dy[2];
00255 gm[0][7]=dx[2];
00256 gm[0][8]=dy[5]+dx[8];
00257 }
00258 }
00259
00260 void planeelemrotlt::addgeommat (matrix &gm,vector &x,vector &y,vector &areacoord,double &jac)
00261 {
00262 long i,i1,i2,i3,ii,jj;
00263 double det,area;
00264 vector b(3),c(3),bf(ndofe),dx(ndofe),dy(ndofe);
00265
00266
00267 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00268 jac = det;
00269 area=det/2.0;
00270 plsb (b.a,y.a,det);
00271 plsc (c.a,x.a,det);
00272
00273 bf_rot_3_2d (bf.a,areacoord.a,b.a,c.a);
00274 dx_bf_rot_3_2d (dx.a,areacoord.a,b.a,c.a);
00275 dy_bf_rot_3_2d (dy.a,areacoord.a,b.a,c.a);
00276
00277 fillm (0.0,gm);
00278
00279 i1=0; i2=1; i3=2; ii=3; jj=6;
00280
00281 for (i=0;i<nne;i++){
00282 gm[0][i1]=0.0-dy[i]/2.0;
00283 gm[0][i2]=dx[i]/2.0;
00284 gm[0][i3]=dx[jj]/2.0-dy[ii]/2.0-bf[i];
00285 i1+=3; i2+=3; i3+=3; ii++; jj++;
00286 }
00287
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 void planeelemrotlt::dmatblock (long ri,long ci,matrix &d, matrix &dd)
00301 {
00302 fillm (0.0,dd);
00303
00304 if (ri==0 && ci==0){
00305 dd[0][0]=d[0][0]; dd[0][1]=d[0][1];
00306 dd[1][0]=d[1][0]; dd[1][1]=d[1][1];
00307 }
00308 if (ri==0 && ci==1){
00309 dd[0][0]=d[0][2];
00310 dd[1][0]=d[1][2];
00311 }
00312 if (ri==1 && ci==0){
00313 dd[0][0]=d[2][0]; dd[0][1]=d[2][1];
00314 }
00315 if (ri==1 && ci==1){
00316 dd[0][0]=d[2][2];
00317 }
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329 void planeelemrotlt::transf_matrix (ivector &nodes,matrix &tmat)
00330 {
00331 long i,n,m;
00332
00333 fillm (0.0,tmat);
00334
00335 n=nodes.n;
00336 m=tmat.m;
00337 for (i=0;i<m;i++){
00338 tmat[i][i]=1.0;
00339 }
00340
00341 for (i=0;i<n;i++){
00342 if (Mt->nodes[nodes[i]].transf>0){
00343 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;
00344 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;
00345 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;
00346 }
00347 }
00348 }
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362 void planeelemrotlt::stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
00363 {
00364 long i,ii,jj,ipp;
00365 double jac,det,thick;
00366 ivector nodes(nne);
00367 vector areacoord(3),t(nne),gp1,gp2,w;
00368 matrix gmr,gmc,dd,d(tncomp,tncomp);
00369
00370 Mt->give_elemnodes (eid,nodes);
00371 Mc->give_thickness (eid,nodes,t);
00372
00373
00374 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00375
00376 fillm (0.0,sm);
00377
00378 for (ii=0;ii<nb;ii++){
00379 allocm (ncomp[ii],ndofe,gmr);
00380 for (jj=0;jj<nb;jj++){
00381 if (intordsm[ii][jj]==0) continue;
00382
00383 allocm (ncomp[jj],ndofe,gmc);
00384 allocm (ncomp[ii],ncomp[jj],dd);
00385
00386 allocv (nip[ii][jj],gp1);
00387 allocv (nip[ii][jj],gp2);
00388 allocv (nip[ii][jj],w);
00389
00390
00391 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
00392 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00393
00394 for (i=0;i<intordsm[ii][jj];i++){
00395 areacoord[0]=gp1[i];
00396 areacoord[1]=gp2[i];
00397 areacoord[2]=1.0-gp1[i]-gp2[i];
00398
00399
00400 geom_matrix_block (gmr,ii,x,y,areacoord);
00401 geom_matrix_block (gmc,jj,x,y,areacoord);
00402
00403
00404
00405
00406
00407
00408
00409 Mm->matstiff (d,ipp);
00410 dmatblock (ii,jj,d,dd);
00411
00412 thick = approx (areacoord,t);
00413
00414 jac=thick*det*w[i];
00415
00416
00417
00418 bdbjac (sm,gmr,dd,gmc,jac);
00419
00420 ipp++;
00421 }
00422
00423 destrv (w);
00424 destrv (gp1);
00425 destrv (gp2);
00426 destrm (gmc);
00427 destrm (dd);
00428 }
00429 destrm (gmr);
00430 }
00431
00432 allocv (1,w);
00433 allocv (1,gp1);
00434 allocv (1,gp2);
00435
00436 gauss_points_tr (gp1.a,gp2.a,w.a,1);
00437 areacoord[0]=gp1[0];
00438 areacoord[1]=gp2[0];
00439 areacoord[2]=1.0-gp1[0]-gp2[0];
00440
00441 allocm (1,ndofe,gmr);
00442 allocm (ndofe,ndofe,gmc);
00443
00444 addgeommat (gmr,x,y,areacoord,jac);
00445 mtxm (gmr,gmr,gmc);
00446 thick = approx (areacoord,t);
00447 cmulm (thick*jac*w[0],gmc);
00448 addm (sm,gmc,sm);
00449
00450 destrm (gmc); destrm (gmr);
00451 destrv (gp1); destrv (gp2); destrv (w);
00452
00453
00454
00455 }
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 void planeelemrotlt::res_stiffness_matrix (long eid,matrix &sm)
00466 {
00467 long transf;
00468 vector x(nne),y(nne);
00469 ivector nodes(nne);
00470 Mt->give_node_coord2d (x,y,eid);
00471 Mt->give_elemnodes (eid,nodes);
00472
00473 stiffness_matrix (eid,0,0,sm,x,y);
00474
00475
00476 transf = Mt->locsystems (nodes);
00477 if (transf>0){
00478 matrix tmat (ndofe,ndofe);
00479 transf_matrix (nodes,tmat);
00480 glmatrixtransf (sm,tmat);
00481 }
00482 }
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495 void planeelemrotlt::mass_matrix (long eid,matrix &mm,vector &x,vector &y)
00496 {
00497 long i;
00498 double jac,area,thick,rho;
00499 ivector nodes(nne);
00500 vector b(3),c(3),areacoord(3),w(intordmm),gp1(intordmm),gp2(intordmm),t(nne),dens(nne);
00501 matrix n(napfun,ndofe);
00502
00503 Mt->give_elemnodes (eid,nodes);
00504 Mt->give_node_coord2d (x,y,eid);
00505 Mc->give_thickness (eid,nodes,t);
00506 Mc->give_density (eid,nodes,dens);
00507
00508
00509 area = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
00510
00511 b_coeff (b.a,y.a);
00512 c_coeff (c.a,x.a);
00513
00514 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00515
00516 fillm (0.0,mm);
00517
00518 for (i=0;i<intordmm;i++){
00519
00520 areacoord[0]=gp1[i];
00521 areacoord[1]=gp2[i];
00522 areacoord[2]=1.0-gp1[i]-gp2[i];
00523
00524 bf_matrix (n,areacoord,b,c);
00525
00526 thick = approx (areacoord,t);
00527 rho = approx (areacoord,dens);
00528
00529 jac=w[i]*thick*rho*area*2.0;
00530
00531 nnj (mm.a,n.a,jac,n.m,n.n);
00532 }
00533
00534 }
00535
00536 void planeelemrotlt::res_mass_matrix (long eid,matrix &mm)
00537 {
00538 vector x(nne),y(nne);
00539 Mt->give_node_coord2d (x,y,eid);
00540 mass_matrix (eid,mm,x,y);
00541 }
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555 void planeelemrotlt::load_matrix (long eid,matrix &lm)
00556 {
00557 long i;
00558 double jac,area,thick;
00559 ivector nodes(nne);
00560 vector x(nne),y(nne),w(intordmm),gp1(intordmm),gp2(intordmm),b(3),c(3),areacoord(3),t(nne);
00561 matrix n(napfun,ndofe);
00562
00563 Mt->give_elemnodes (eid,nodes);
00564 Mt->give_node_coord2d (x,y,eid);
00565 Mc->give_thickness (eid,nodes,t);
00566
00567 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00568
00569 area = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
00570 b_coeff (b.a,y.a);
00571 c_coeff (c.a,x.a);
00572
00573 fillm (0.0,lm);
00574
00575 for (i=0;i<intordmm;i++){
00576 areacoord[0]=gp1[i];
00577 areacoord[1]=gp2[i];
00578 areacoord[2]=1.0-gp1[i]-gp2[i];
00579
00580 bf_matrix (n,areacoord,b,c);
00581
00582 thick = approx (areacoord,t);
00583
00584
00585 jac=w[i]*thick*area*2.0;
00586
00587 nnj (lm.a,n.a,jac,n.m,n.n);
00588 }
00589
00590 }
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 void planeelemrotlt::res_ip_strains (long lcid,long eid)
00601 {
00602 vector aux,x(nne),y(nne),r(ndofe);
00603 ivector nodes(nne);
00604 matrix tmat;
00605
00606 Mt->give_node_coord2d (x,y,eid);
00607 Mt->give_elemnodes (eid,nodes);
00608 eldispl (lcid,eid,r.a);
00609
00610
00611 long transf = Mt->locsystems (nodes);
00612 if (transf>0){
00613 allocv (ndofe,aux);
00614 allocm (ndofe,ndofe,tmat);
00615 transf_matrix (nodes,tmat);
00616
00617 lgvectortransf (aux,r,tmat);
00618 copyv (aux,r);
00619 destrv (aux);
00620 destrm (tmat);
00621 }
00622
00623 ip_strains (lcid,eid,0,0,x,y,r);
00624
00625 }
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638 void planeelemrotlt::ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00639 {
00640 long i,ii,ipp;
00641 vector gp1,gp2,w,eps,epscum,areacoord(3);
00642 matrix gm;
00643
00644 allocm (tncomp,ndofe,gm);
00645 allocv (tncomp,eps);
00646
00647
00648 for (ii=0;ii<nb;ii++){
00649
00650 allocv (intordsm[ii][ii],gp1);
00651 allocv (intordsm[ii][ii],gp2);
00652 allocv (intordsm[ii][ii],w);
00653
00654
00655
00656 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00657
00658 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00659 for (i=0;i<intordsm[ii][ii];i++){
00660 areacoord[0]=gp1[i];
00661 areacoord[1]=gp2[i];
00662 areacoord[2]=1.0-gp1[i]-gp2[i];
00663
00664
00665
00666 geom_matrix (gm,x,y,areacoord);
00667
00668 mxv (gm,r,eps);
00669
00670
00671 Mm->storestrain (lcid,ipp,0,tncomp,eps);
00672 ipp++;
00673 }
00674
00675
00676 destrv (w); destrv (gp1); destrv (gp2);
00677 }
00678 destrm (gm); destrv (eps);
00679 }
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690 void planeelemrotlt::nod_strains_ip (long lcid,long eid,long ri,long ci)
00691 {
00692 long i,j,ipp;
00693 ivector ipnum(nne),nod(nne);
00694 vector eps(gncomp);
00695
00696
00697
00698 ipp=Mt->elements[eid].ipp[ri][ci];
00699 nodip_planelt (ipp,intordsm[0][0],ipnum);
00700
00701
00702 Mt->give_elemnodes (eid,nod);
00703
00704 for (i=0;i<nne;i++){
00705
00706 Mm->givestrain (lcid,ipnum[i],eps);
00707
00708
00709 j=nod[i];
00710 Mt->nodes[j].storestrain (lcid,0,eps);
00711 }
00712
00713 }
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725 void planeelemrotlt::ip_stresses (long lcid,long eid,long ri,long ci)
00726 {
00727 compute_nlstress (lcid,eid,ri,ci);
00728 }
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738 void planeelemrotlt::res_ip_stresses (long lcid,long eid)
00739 {
00740 long ri,ci;
00741 ri=0;
00742 ci=0;
00743 compute_nlstress (lcid,eid,ri,ci);
00744 }
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755 void planeelemrotlt::nod_stresses_ip (long lcid,long eid,long ri,long ci)
00756 {
00757 long i,j,ipp;
00758 ivector ipnum(nne),nod(nne);
00759 vector sig(gncomp);
00760
00761
00762
00763 ipp=Mt->elements[eid].ipp[ri][ci];
00764 nodip_planelt (ipp,intordsm[0][0],ipnum);
00765
00766
00767 Mt->give_elemnodes (eid,nod);
00768
00769 for (i=0;i<nne;i++){
00770
00771 Mm->givestress (lcid,ipnum[i],sig);
00772
00773
00774 j=nod[i];
00775 Mt->nodes[j].storestress (lcid,0,sig);
00776 }
00777
00778 }
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789 void planeelemrotlt::compute_nlstress (long lcid,long eid,long ri,long ci)
00790 {
00791 long i,ii,ipp;
00792
00793 for (ii=0;ii<nb;ii++){
00794 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00795
00796 for (i=0;i<intordsm[ii][ii];i++){
00797
00798 if (Mp->strcomp==1)
00799 Mm->computenlstresses (ipp);
00800
00801 ipp++;
00802 }
00803 }
00804 }
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815 void planeelemrotlt::compute_nlstressincr (long lcid,long eid,long ri,long ci)
00816 {
00817 long i,ii,ipp;
00818
00819 for (ii=0;ii<nb;ii++){
00820 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00821
00822 for (i=0;i<intordsm[ii][ii];i++){
00823
00824 if (Mp->strcomp==1)
00825 Mm->computenlstressesincr (ipp);
00826
00827 ipp++;
00828 }
00829 }
00830 }
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847 void planeelemrotlt::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00848 {
00849 integratedquant iq;
00850 iq=locstress;
00851
00852
00853 compute_nlstress (lcid,eid,ri,ci);
00854
00855
00856 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00857
00858 }
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883 void planeelemrotlt::res_mainip_strains (long lcid,long eid)
00884 {
00885 vector aux,x(nne),y(nne),r(ndofe);
00886 ivector nodes(nne);
00887 matrix tmat;
00888
00889 Mt->give_node_coord2d (x,y,eid);
00890 Mt->give_elemnodes (eid,nodes);
00891 eldispl (lcid,eid,r.a);
00892
00893
00894 long transf = Mt->locsystems (nodes);
00895 if (transf>0){
00896 allocv (ndofe,aux);
00897 allocm (ndofe,ndofe,tmat);
00898 transf_matrix (nodes,tmat);
00899
00900 lgvectortransf (aux,r,tmat);
00901 copyv (aux,r);
00902 destrv (aux);
00903 destrm (tmat);
00904 }
00905
00906 mainip_strains (lcid,eid,0,0,x,y,r);
00907 }
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920 void planeelemrotlt::mainip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00921 {
00922 long i,ii,ipp;
00923 vector gp1,gp2,w,eps,areacoord(3);
00924 matrix gm;
00925
00926 for (ii=0;ii<nb;ii++){
00927 allocv (intordsm[ii][ii],gp1);
00928 allocv (intordsm[ii][ii],gp2);
00929 allocv (intordsm[ii][ii],w);
00930 allocv (ncomp[ii],eps);
00931 allocm (ncomp[ii],ndofe,gm);
00932
00933 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00934
00935 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00936 for (i=0;i<intordsm[ii][ii];i++){
00937 areacoord[0]=gp1[i];
00938 areacoord[1]=gp2[i];
00939 areacoord[2]=1.0-gp1[i]-gp2[i];
00940
00941 geom_matrix_block (gm,ii,x,y,areacoord);
00942 mxv (gm,r,eps);
00943
00944 Mm->storestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00945 ipp++;
00946 }
00947
00948 destrm (gm); destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
00949 }
00950 }
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960 void planeelemrotlt::nod_strains (long lcid,long eid,long ri,long ci)
00961 {
00962 long i,ii,ipp;
00963 double *lsm,*lhs,*rhs;
00964 vector nxi(nne),neta(nne),gp1,gp2,w,eps,aux,natcoord(2),areacoord(3);
00965 ivector nodes(nne);
00966
00967 lsm = new double [9];
00968
00969 nodecoord (nxi,neta);
00970 Mt->give_elemnodes (eid,nodes);
00971
00972
00973 for (ii=0;ii<nb;ii++){
00974 allocv (intordsm[ii][ii],gp1);
00975 allocv (intordsm[ii][ii],gp2);
00976 allocv (intordsm[ii][ii],w);
00977 allocv (ncomp[ii],eps);
00978 lhs = new double [ncomp[ii]*3];
00979 rhs = new double [ncomp[ii]*3];
00980 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00981
00982 nullv (lsm,9);
00983 nullv (rhs,ncomp[ii]*3);
00984
00985 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00986 for (i=0;i<intordsm[ii][ii];i++){
00987 areacoord[0]=gp1[i];
00988 areacoord[1]=gp2[i];
00989 areacoord[2]=1.0-gp1[i]-gp2[i];
00990
00991 Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00992
00993 natcoord[0]=areacoord[0]; natcoord[1]=areacoord[1];
00994 matassem_lsm (lsm,natcoord);
00995 rhsassem_lsm (rhs,natcoord,eps);
00996
00997 ipp++;
00998 }
00999
01000 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
01001 Mt->strain_nodal_values (nodes,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii],lcid);
01002
01003 delete [] lhs; delete [] rhs;
01004 destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
01005 }
01006
01007 delete [] lsm;
01008 }
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019 void planeelemrotlt::elem_strains (double **stra,long lcid,long eid,long ri,long ci)
01020 {
01021 long i,ii,ipp;
01022 double *lsm,*lhs,*rhs;
01023 vector nxi(nne),neta(nne),gp1,gp2,w,eps,aux,natcoord(2),areacoord(3);
01024 ivector nodes(nne);
01025
01026 lsm = new double [9];
01027 nodecoord (nxi,neta);
01028 Mt->give_elemnodes (eid,nodes);
01029
01030 for (ii=0;ii<nb;ii++){
01031 allocv (intordsm[ii][ii],gp1);
01032 allocv (intordsm[ii][ii],gp2);
01033 allocv (intordsm[ii][ii],w);
01034 allocv (ncomp[ii],eps);
01035 lhs = new double [ncomp[ii]*3];
01036 rhs = new double [ncomp[ii]*3];
01037 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01038
01039 nullv (lsm,9);
01040 nullv (rhs,ncomp[ii]*3);
01041
01042 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01043 for (i=0;i<intordsm[ii][ii];i++){
01044 areacoord[0]=gp1[i];
01045 areacoord[1]=gp2[i];
01046 areacoord[2]=1.0-areacoord[0]-areacoord[1];
01047
01048 Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
01049
01050 natcoord[0]=areacoord[0]; natcoord[1]=areacoord[1];
01051 matassem_lsm (lsm,natcoord);
01052 rhsassem_lsm (rhs,natcoord,eps);
01053
01054 ipp++;
01055 }
01056
01057 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
01058 nodal_values (stra,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii]);
01059
01060 delete [] lhs; delete [] rhs;
01061 destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
01062 }
01063
01064 delete [] lsm;
01065 }
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077 void planeelemrotlt::appstrain (long lcid,long eid,double xi,double eta,long fi,long ncomp,vector &eps)
01078 {
01079 long i,j,k;
01080 ivector nodes(nne);
01081 vector nodval(nne),areacoord(3);
01082
01083 if (ncomp != eps.n){
01084 fprintf (stderr,"\n\n wrong interval of indices in function strain (%s, line %d).\n",__FILE__,__LINE__);
01085 abort ();
01086 }
01087
01088 areacoord[0]=xi;
01089 areacoord[1]=eta;
01090 areacoord[2]=1.0-areacoord[0]-areacoord[1];
01091
01092 Mt->give_elemnodes (eid,nodes);
01093 k=0;
01094 for (i=fi;i<fi+ncomp;i++){
01095 for (j=0;j<nne;j++){
01096 nodval[j]=Mt->nodes[nodes[j]].strain[lcid*tncomp+i];
01097 }
01098 eps[k]=approx (areacoord,nodval);
01099 k++;
01100 }
01101
01102 }
01103
01104 void planeelemrotlt::allip_strains (double **stra,long lcid,long eid,long ri,long ci)
01105 {
01106 long i,ii,jj,ipp;
01107 vector eps(tncomp),gp1,gp2,w;
01108
01109 allocv (tncomp,eps);
01110
01111 for (ii=0;ii<nb;ii++){
01112 for (jj=0;jj<nb;jj++){
01113 if (intordsm[ii][jj]==0) continue;
01114
01115 allocv (intordsm[ii][jj],gp1);
01116 allocv (intordsm[ii][jj],gp2);
01117 allocv (intordsm[ii][jj],w);
01118 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
01119 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01120
01121 for (i=0;i<intordsm[ii][jj];i++){
01122
01123 if (Mp->strainaver==0)
01124 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
01125 if (Mp->strainaver==1)
01126 appstrain (lcid,eid,gp1[i],gp2[i],0,tncomp,eps);
01127
01128 Mm->storestrain (lcid,ipp,eps);
01129 ipp++;
01130 }
01131 destrv (w); destrv (gp2); destrv (gp1);
01132 }
01133 }
01134 destrv (eps);
01135 }
01136
01137 void planeelemrotlt::strains (long lcid,long eid,long ri,long ci)
01138 {
01139 long i,naep,ncp,sid;
01140 double **stra;
01141 vector coord,eps;
01142
01143 if (Mp->strainaver==0){
01144 stra = new double* [nne];
01145 for (i=0;i<nne;i++){
01146 stra[i] = new double [tncomp];
01147 }
01148 elem_strains (stra,lcid,eid,ri,ci);
01149 }
01150
01151 switch (Mm->stra.tape[eid]){
01152 case nowhere:{
01153 break;
01154 }
01155 case intpts:{
01156 allip_strains (stra,lcid,eid,ri,ci);
01157 break;
01158 }
01159 case enodes:{
01160 break;
01161 }
01162 case userdefined:{
01163
01164 naep = Mm->stra.give_naep (eid);
01165 ncp = Mm->stra.give_ncomp (eid);
01166 sid = Mm->stra.give_sid (eid);
01167 allocv (ncp,eps);
01168 allocv (2,coord);
01169 for (i=0;i<naep;i++){
01170 Mm->stra.give_aepcoord (sid,i,coord);
01171
01172 if (Mp->strainaver==0)
01173 appval (coord[0],coord[1],0,ncp,eps,stra);
01174 if (Mp->strainaver==1)
01175 appstrain (lcid,eid,coord[0],coord[1],0,ncp,eps);
01176
01177 Mm->stra.storevalues(lcid,eid,i,eps);
01178 }
01179 destrv (eps);
01180 destrv (coord);
01181 break;
01182 }
01183 default:{
01184 fprintf (stderr,"\n\n unknown strain point is required in function planeelemrotlt::strains (%s, line %d).\n",__FILE__,__LINE__);
01185 }
01186 }
01187
01188 if (Mp->strainaver==0){
01189 for (i=0;i<nne;i++){
01190 delete [] stra[i];
01191 }
01192 delete [] stra;
01193 }
01194 }
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204 void planeelemrotlt::nodecoord (vector &xi,vector &eta)
01205 {
01206 xi[0] = 0.0; eta[0] = 0.0;
01207 xi[1] = 1.0; eta[1] = 0.0;
01208 xi[2] = 0.0; eta[2] = 1.0;
01209 }
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220 void planeelemrotlt::appval (double xi,double eta,long fi,long nc,vector &eps,double **val)
01221 {
01222 long i,j,k;
01223 vector nodval;
01224
01225 k=0;
01226 allocv (nne,nodval);
01227 for (i=fi;i<fi+nc;i++){
01228 for (j=0;j<nne;j++){
01229 nodval[j]=val[j][i];
01230 }
01231 eps[k]=approx_nat (xi,eta,nodval);
01232 k++;
01233 }
01234
01235 destrv (nodval);
01236 }
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247 void planeelemrotlt::mainip_stresses (long lcid,long eid,long ri,long ci)
01248 {
01249 long i,ii,jj,ipp;
01250 vector gp1,gp2,w,eps,sig,auxsig;
01251 matrix d(tncomp,tncomp),dd;
01252
01253 for (ii=0;ii<nb;ii++){
01254 allocv (ncomp[ii],sig);
01255 allocv (ncomp[ii],auxsig);
01256 allocv (intordsm[ii][ii],gp1);
01257 allocv (intordsm[ii][ii],gp2);
01258 allocv (intordsm[ii][ii],w);
01259
01260 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01261 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01262
01263 for (i=0;i<intordsm[ii][ii];i++){
01264
01265 Mm->matstiff (d,ipp);
01266
01267 fillv (0.0,sig);
01268 for (jj=0;jj<nb;jj++){
01269 allocv (ncomp[jj],eps);
01270 allocm (ncomp[ii],ncomp[jj],dd);
01271
01272 if (Mp->strainaver==0)
01273 Mm->givestrain (lcid,ipp,cncomp[jj],ncomp[jj],eps);
01274 if (Mp->strainaver==1)
01275 appstrain (lcid,eid,gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps);
01276
01277 dmatblock (ii,jj,d,dd);
01278 mxv (dd,eps,auxsig);
01279 addv (auxsig,sig,sig);
01280 destrm (dd); destrv (eps);
01281 }
01282
01283 Mm->storestress (lcid,ipp,cncomp[ii],ncomp[ii],sig);
01284
01285 ipp++;
01286 }
01287 destrv (w); destrv (gp1); destrv (gp2); destrv (auxsig); destrv (sig);
01288 }
01289
01290 }
01291
01292 void planeelemrotlt::nod_stresses (long lcid,long eid,long ri,long ci)
01293 {
01294 long i,ii,jj,ipp;
01295 double *lsm,*lhs,*rhs;
01296 vector nxi(nne),neta(nne),gp1,gp2,w,eps,sig,auxsig,natcoord(2);
01297 ivector nodes(nne);
01298 matrix d(tncomp,tncomp),dd;
01299
01300 lsm = new double [9];
01301
01302 nodecoord (nxi,neta);
01303 Mt->give_elemnodes (eid,nodes);
01304
01305 for (ii=0;ii<nb;ii++){
01306 allocv (intordsm[ii][ii],gp1);
01307 allocv (intordsm[ii][ii],gp2);
01308 allocv (intordsm[ii][ii],w);
01309 allocv (ncomp[ii],sig);
01310 allocv (ncomp[ii],auxsig);
01311 lhs = new double [ncomp[ii]*3];
01312 rhs = new double [ncomp[ii]*3];
01313 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01314
01315 nullv (lsm,9);
01316 nullv (rhs,ncomp[ii]*3);
01317
01318 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01319
01320 for (i=0;i<intordsm[ii][ii];i++){
01321
01322 Mm->matstiff (d,ipp);
01323 ipp++;
01324
01325 fillv (0.0,sig);
01326 for (jj=0;jj<nb;jj++){
01327 allocv (ncomp[jj],eps);
01328 allocm (ncomp[ii],ncomp[jj],dd);
01329
01330 if (Mp->strainaver==0)
01331 Mm->givestrain (lcid,ipp,cncomp[jj],ncomp[jj],eps);
01332 if (Mp->strainaver==1)
01333 appstrain (lcid,eid,gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps);
01334
01335 dmatblock (ii,jj,d,dd);
01336 mxv (dd,eps,auxsig);
01337 addv (auxsig,sig,sig);
01338 destrm (dd); destrv (eps);
01339 }
01340
01341 natcoord[0]=gp1[i]; natcoord[1]=gp2[i];
01342 matassem_lsm (lsm,natcoord);
01343 rhsassem_lsm (rhs,natcoord,sig);
01344
01345 }
01346
01347 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
01348 Mt->stress_nodal_values (nodes,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii],lcid);
01349
01350
01351 delete [] lhs; delete [] rhs;
01352 destrv (auxsig); destrv (sig); destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
01353 }
01354
01355 delete [] lsm;
01356 }
01357
01358 void planeelemrotlt::elem_stresses (double **stra,double **stre,long lcid,long eid,long ri,long ci)
01359 {
01360 long i,ii,jj,ipp;
01361 double *lsm,*lhs,*rhs;
01362 vector gp1,gp2,w,nxi(nne),neta(nne),eps,sig,auxsig,natcoord(2);
01363 ivector nodes(nne);
01364 matrix d(tncomp,tncomp),dd;
01365
01366 lsm = new double [9];
01367
01368 nodecoord (nxi,neta);
01369 Mt->give_elemnodes (eid,nodes);
01370
01371 for (ii=0;ii<nb;ii++){
01372 allocv (intordsm[ii][ii],gp1);
01373 allocv (intordsm[ii][ii],gp2);
01374 allocv (intordsm[ii][ii],w);
01375 allocv (ncomp[ii],sig);
01376 allocv (ncomp[ii],auxsig);
01377 lhs = new double [ncomp[ii]*3];
01378 rhs = new double [ncomp[ii]*3];
01379 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01380
01381 nullv (lsm,9);
01382 nullv (rhs,ncomp[ii]*3);
01383
01384 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01385
01386 for (i=0;i<intordsm[ii][ii];i++){
01387 Mm->matstiff (d,ipp);
01388 ipp++;
01389
01390 fillv (0.0,sig);
01391 for (jj=0;jj<nb;jj++){
01392 allocv (ncomp[jj],eps);
01393 allocm (ncomp[ii],ncomp[jj],dd);
01394
01395 if (Mp->strainaver==0)
01396 appval (gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps,stra);
01397 if (Mp->strainaver==1)
01398 appstrain (lcid,eid,gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps);
01399
01400 dmatblock (ii,jj,d,dd);
01401 mxv (dd,eps,auxsig);
01402 addv (auxsig,sig,sig);
01403 destrm (dd); destrv (eps);
01404 }
01405
01406 natcoord[0]=gp1[i]; natcoord[1]=gp2[i];
01407 matassem_lsm (lsm,natcoord);
01408 rhsassem_lsm (rhs,natcoord,sig);
01409 }
01410
01411 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
01412 nodal_values (stre,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii]);
01413
01414
01415 delete [] lhs; delete [] rhs;
01416 destrv (sig); destrv (eps);
01417 }
01418
01419 delete [] lsm;
01420 }
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432 void planeelemrotlt::appstress (long lcid,long eid,double xi,double eta,long fi,long ncomp,vector &sig)
01433 {
01434 long i,j,k;
01435 ivector nodes(nne);
01436 vector areacoord(3),nodval(nne);
01437
01438 if (ncomp != sig.n){
01439 fprintf (stderr,"\n\n wrong interval of indices in function stress (%s, line %d).\n",__FILE__,__LINE__);
01440 abort ();
01441 }
01442
01443 areacoord[0]=xi;
01444 areacoord[1]=eta;
01445 areacoord[2]=1.0-areacoord[0]-areacoord[1];
01446
01447 Mt->give_elemnodes (eid,nodes);
01448 k=0;
01449 for (i=fi;i<fi+ncomp;i++){
01450 for (j=0;j<nne;j++){
01451 nodval[j]=Mt->nodes[nodes[j]].stress[lcid*tncomp+i];
01452 }
01453 sig[k]=approx (areacoord,nodval);
01454 k++;
01455 }
01456
01457 }
01458
01459 void planeelemrotlt::allip_stresses (double **stre,long lcid,long eid,long ri,long ci)
01460 {
01461 long i,ii,jj,ipp;
01462 vector sig(tncomp),gp1,gp2,w;
01463
01464 for (ii=0;ii<nb;ii++){
01465 for (jj=0;jj<nb;jj++){
01466 if (intordsm[ii][jj]==0) continue;
01467 allocv (intordsm[ii][jj],gp1);
01468 allocv (intordsm[ii][jj],gp2);
01469 allocv (intordsm[ii][jj],w);
01470 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
01471 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01472
01473 for (i=0;i<intordsm[ii][jj];i++){
01474
01475 if (Mp->stressaver==0)
01476 appval (gp1[i],gp2[i],0,tncomp,sig,stre);
01477 if (Mp->stressaver==1)
01478 appstress (lcid,eid,gp1[i],gp2[i],0,tncomp,sig);
01479
01480 Mm->storestress (lcid,ipp,sig);
01481 ipp++;
01482 }
01483 }
01484 destrv (w); destrv (gp1); destrv (gp2);
01485 }
01486 }
01487
01488 void planeelemrotlt::stresses (long lcid,long eid,long ri,long ci)
01489 {
01490 long i,naep,ncp,sid;
01491 double **stra,**stre;
01492 vector coord,sig;
01493
01494 if (Mp->stressaver==0){
01495 stra = new double* [nne];
01496 stre = new double* [nne];
01497 for (i=0;i<nne;i++){
01498 stra[i] = new double [tncomp];
01499 stre[i] = new double [tncomp];
01500 }
01501 elem_strains (stra,lcid,eid,ri,ci);
01502 elem_stresses (stra,stre,lcid,eid,ri,ci);
01503 }
01504
01505 switch (Mm->stre.tape[eid]){
01506 case nowhere:{
01507 break;
01508 }
01509 case intpts:{
01510 allip_stresses (stre,lcid,eid,ri,ci);
01511 break;
01512 }
01513 case enodes:{
01514 break;
01515 }
01516 case userdefined:{
01517
01518 naep = Mm->stre.give_naep (eid);
01519 ncp = Mm->stre.give_ncomp (eid);
01520 sid = Mm->stre.give_sid (eid);
01521 allocv (ncp,sig);
01522 allocv (2,coord);
01523 for (i=0;i<naep;i++){
01524 Mm->stre.give_aepcoord (sid,i,coord);
01525
01526 if (Mp->stressaver==0)
01527 appval (coord[0],coord[1],0,ncp,sig,stre);
01528 if (Mp->stressaver==1)
01529 appstress (lcid,eid,coord[0],coord[1],0,ncp,sig);
01530
01531 Mm->stre.storevalues(lcid,eid,i,sig);
01532 }
01533 destrv (sig);
01534 destrv (coord);
01535 break;
01536 }
01537 default:{
01538 fprintf (stderr,"\n\n unknown stress point is required in function planeelemrotlt::stresses (%s, line %d).\n",__FILE__,__LINE__);
01539 }
01540 }
01541
01542 if (Mp->stressaver==0){
01543 for (i=0;i<nne;i++){
01544 delete [] stra[i];
01545 delete [] stre[i];
01546 }
01547 delete [] stra;
01548 delete [] stre;
01549 }
01550 }
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561 void planeelemrotlt::nod_others (long lcid,long eid,long ri,long ci)
01562 {
01563 long i,ii,j,ncomp,ipp;
01564 double *lsm,*lhs,*rhs;
01565 vector nxi(nne),neta(nne),gp1,gp2,w,other,aux,natcoord(2),areacoord(3);
01566 ivector nodes(nne);
01567
01568 lsm = new double [9];
01569
01570 nodecoord (nxi,neta);
01571 Mt->give_elemnodes (eid,nodes);
01572
01573
01574 for (ii=0;ii<nb;ii++){
01575 allocv (intordsm[ii][ii],gp1);
01576 allocv (intordsm[ii][ii],gp2);
01577 allocv (intordsm[ii][ii],w);
01578 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01579 ncomp = Mm->ip[ipp].ncompother;
01580 allocv (ncomp,other);
01581 lhs = new double [ncomp*3];
01582 rhs = new double [ncomp*3];
01583 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01584
01585 nullv (lsm,9);
01586 nullv (rhs,ncomp*3);
01587
01588 for (i=0;i<intordsm[ii][ii];i++){
01589 areacoord[0]=gp1[i];
01590 areacoord[1]=gp2[i];
01591 areacoord[2]=1.0-gp1[i]-gp2[i];
01592
01593 for (j = 0; j < ncomp; j++)
01594 other[j] = Mm->ip[ipp].eqother[j];
01595
01596 natcoord[0]=areacoord[0]; natcoord[1]=areacoord[1];
01597 matassem_lsm (lsm,natcoord);
01598 rhsassem_lsm (rhs,natcoord,other);
01599
01600 ipp++;
01601 }
01602
01603 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp);
01604 Mt->other_nodal_values (nodes,nxi,neta,nxi,lhs,2,0,ncomp);
01605
01606 delete [] lhs; delete [] rhs;
01607 destrv (other); destrv (w); destrv (gp1); destrv (gp2);
01608 }
01609
01610 delete [] lsm;
01611 }
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630 void planeelemrotlt::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01631 {
01632 integratedquant iq;
01633 iq=nonlocstress;
01634
01635
01636 compute_nonloc_nlstress (lcid,eid,ri,ci);
01637
01638 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01639 }
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654 void planeelemrotlt::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01655 {
01656 integratedquant iq;
01657 iq=stressincr;
01658
01659
01660 compute_nlstressincr (lcid,eid,ri,ci);
01661
01662 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01663 }
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678 void planeelemrotlt::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y)
01679 {
01680 integratedquant iq;
01681 iq=eigstress;
01682
01683
01684 if (Mp->eigstrcomp)
01685 compute_eigstress (lcid,eid,ri,ci);
01686
01687 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y);
01688 }
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701 void planeelemrotlt::res_internal_forces (long lcid,long eid,vector &ifor)
01702 {
01703 long transf;
01704 ivector nodes (nne);
01705 vector v(ndofe),x(nne),y(nne);
01706
01707 Mt->give_node_coord2d (x,y,eid);
01708 internal_forces (lcid,eid,0,0,ifor,x,y);
01709
01710
01711 Mt->give_elemnodes (eid,nodes);
01712 transf = Mt->locsystems (nodes);
01713 if (transf>0){
01714 matrix tmat (ndofe,ndofe);
01715 transf_matrix (nodes,tmat);
01716 glvectortransf (ifor,v,tmat);
01717 copyv (v,ifor);
01718 }
01719 }
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732 void planeelemrotlt::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
01733 {
01734 long transf;
01735 ivector nodes (nne);
01736 vector v(ndofe),x(nne),y(nne);
01737
01738 Mt->give_node_coord2d (x,y,eid);
01739 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y);
01740
01741
01742 Mt->give_elemnodes (eid,nodes);
01743 transf = Mt->locsystems (nodes);
01744 if (transf>0){
01745 matrix tmat (ndofe,ndofe);
01746 transf_matrix (nodes,tmat);
01747 glvectortransf (ifor,v,tmat);
01748 copyv (v,ifor);
01749 }
01750 }
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763 void planeelemrotlt::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01764 {
01765 long transf;
01766 ivector nodes (nne);
01767 vector v(ndofe),x(nne),y(nne);
01768
01769 Mt->give_node_coord2d (x,y,eid);
01770 incr_internal_forces (lcid,eid,0,0,ifor,x,y);
01771
01772
01773 Mt->give_elemnodes (eid,nodes);
01774 transf = Mt->locsystems (nodes);
01775 if (transf>0){
01776 matrix tmat (ndofe,ndofe);
01777 transf_matrix (nodes,tmat);
01778 glvectortransf (ifor,v,tmat);
01779 copyv (v,ifor);
01780 }
01781 }
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794 void planeelemrotlt::res_eigstrain_forces (long lcid,long eid,vector &nfor)
01795 {
01796 long transf;
01797 ivector nodes (nne);
01798 vector v(ndofe),x(nne),y(nne);
01799
01800 Mt->give_node_coord2d (x,y,eid);
01801 eigstrain_forces (lcid,eid,0,0,nfor,x,y);
01802 Mt->give_node_coord2d (x,y,eid);
01803 internal_forces (lcid,eid,0,0,nfor,x,y);
01804
01805
01806 Mt->give_elemnodes (eid,nodes);
01807 transf = Mt->locsystems (nodes);
01808 if (transf>0){
01809 matrix tmat (ndofe,ndofe);
01810 transf_matrix (nodes,tmat);
01811 glvectortransf (nfor,v,tmat);
01812 copyv (v,nfor);
01813 }
01814 }
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831 void planeelemrotlt::local_values (long lcid,long eid,long ri,long ci)
01832 {
01833 long i,ii,ipp;
01834 double **stra;
01835 ivector nodes(nne);
01836 vector w,gp1,gp2,eps(tncomp);
01837
01838 stra = new double* [nne];
01839 for (i=0;i<nne;i++){
01840 stra[i] = new double [tncomp];
01841 }
01842 elem_strains (stra,lcid,eid,ri,ci);
01843
01844 for (ii=0;ii<nb;ii++){
01845 allocv (intordsm[ii][ii],gp1);
01846 allocv (intordsm[ii][ii],gp2);
01847 allocv (intordsm[ii][ii],w);
01848 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01849 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01850
01851
01852 for (i=0;i<intordsm[ii][ii];i++){
01853 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
01854 Mm->storestrain (lcid,ipp,eps);
01855 if (Mp->strcomp==1)
01856 Mm->computenlstresses (ipp);
01857 ipp++;
01858 }
01859 destrv (w); destrv (gp1); destrv (gp2);
01860 }
01861 for (i=0;i<nne;i++){
01862 delete [] stra[i];
01863 }
01864 delete [] stra;
01865 }
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878 void planeelemrotlt::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01879 {
01880 long i,ii,ipp;
01881 double **stra;
01882 ivector nodes(nne);
01883 vector w,gp1,gp2,eps(tncomp);
01884
01885 stra = new double* [nne];
01886 for (i=0;i<nne;i++){
01887 stra[i] = new double [tncomp];
01888 }
01889 elem_strains (stra,lcid,eid,ri,ci);
01890
01891 for (ii=0;ii<nb;ii++){
01892 allocv (intordsm[ii][ii],gp1);
01893 allocv (intordsm[ii][ii],gp2);
01894 allocv (intordsm[ii][ii],w);
01895 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01896 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01897
01898
01899 for (i=0;i<intordsm[ii][ii];i++){
01900 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
01901 Mm->storestrain (lcid,ipp,eps);
01902 if (Mp->strcomp==1)
01903 Mm->compnonloc_nlstresses (ipp);
01904 ipp++;
01905 }
01906 destrv (w); destrv (gp1); destrv (gp2);
01907 }
01908 for (i=0;i<nne;i++){
01909 delete [] stra[i];
01910 }
01911 delete [] stra;
01912 }
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925 void planeelemrotlt::compute_eigstress (long lcid,long eid,long ri,long ci)
01926 {
01927 long i,ii,ipp;
01928 vector eigstr(tncomp),sig(tncomp);
01929 matrix d(tncomp,tncomp);
01930
01931 for (ii=0;ii<nb;ii++){
01932 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01933
01934 for (i=0;i<intordsm[ii][ii];i++){
01935 Mm->giveeigstrain (ipp,eigstr);
01936
01937 Mm->matstiff (d,ipp);
01938 mxv (d,eigstr,sig);
01939 Mm->storeeigstress (ipp,sig);
01940 ipp++;
01941 }
01942 }
01943 }
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960 void planeelemrotlt::elem_integration (integratedquant iq, long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y)
01961 {
01962 long i,ii,ipp;
01963 double det,thick;
01964 ivector nodes(nne);
01965 vector areacoord(3),t(nne),ipv,contr(ndofe);
01966 vector w,gp1,gp2;
01967 matrix gm;
01968
01969 Mt->give_elemnodes (eid,nodes);
01970 Mc->give_thickness (eid,nodes,t);
01971
01972 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]));
01973
01974 fillv (0.0,nv);
01975
01976 for (ii=0;ii<nb;ii++){
01977 allocv (intordsm[ii][ii],gp1);
01978 allocv (intordsm[ii][ii],gp2);
01979 allocv (intordsm[ii][ii],w);
01980 allocm (ncomp[ii],ndofe,gm);
01981 allocv (ncomp[ii],ipv);
01982
01983 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01984
01985 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01986
01987
01988 for (i=0;i<intordsm[ii][ii];i++){
01989 areacoord[0]=gp1[i];
01990 areacoord[1]=gp2[i];
01991 areacoord[2]=1.0-gp1[i]-gp2[i];
01992
01993 thick = approx (areacoord,t);
01994
01995
01996 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
01997
01998
01999 geom_matrix_block (gm,ii,x,y,areacoord);
02000
02001 mtxv (gm,ipv,contr);
02002
02003 cmulv (det*w[i]*thick,contr);
02004
02005
02006 addv(contr,nv,nv);
02007 ipp++;
02008 }
02009 destrv (w); destrv (gp1); destrv (gp2);
02010 destrm (gm); destrv (ipv);
02011 }
02012 }
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026 void planeelemrotlt::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
02027 {
02028 long i,ii;
02029 vector x(nne),y(nne),areacoord(3),w(intordsm[ri][ci]),gp1(intordsm[ri][ci]),gp2(intordsm[ri][ci]);
02030
02031 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ri][ci]);
02032 Mt->give_node_coord2d (x,y,eid);
02033 ii=Mt->elements[eid].ipp[ri][ci];
02034
02035 for (i=0;i<intordsm[ri][ci];i++){
02036 areacoord[0]=gp1[i];
02037 areacoord[1]=gp2[i];
02038 areacoord[2]=1.0-areacoord[0]-areacoord[1];
02039
02040 if (ii==ipp){
02041 coord[0]=approx (areacoord,x);
02042 coord[1]=approx (areacoord,y);
02043 coord[2]=0.0;
02044 }
02045 ii++;
02046 }
02047 }
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073 void planeelemrotlt::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
02074 {
02075 long i, j, k, ipp;
02076 long ii, jj, nv = nodval.n;
02077 double xi, eta, ipval;
02078 vector w, gp1, gp2, anv(nne);
02079 long nstra, nstre, ncompstr, ncompeqother;
02080 long idstra, idstre, idoth, idic;
02081 inictype ict;
02082
02083 nstra = idstra = nstre = idstre = idoth = idic = 0;
02084
02085 ict = ictn[0];
02086 for (i=0; i<nne; i++)
02087 {
02088 if (ictn[i] != ict)
02089 {
02090 print_err("Incompatible types of initial conditions on element %ld\n"
02091 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
02092 abort();
02093 }
02094 }
02095 for (j = 0; j < nv; j++)
02096 {
02097 for(i = 0; i < nne; i++)
02098 anv[i] = nodval[i][j];
02099 for (ii = 0; ii < nb; ii++)
02100 {
02101 for (jj = 0; jj < nb; jj++)
02102 {
02103 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
02104 if (intordsm[ii][jj] == 0)
02105 continue;
02106 allocv (intordsm[ii][jj],gp1);
02107 allocv (intordsm[ii][jj],gp2);
02108 allocv (intordsm[ii][jj],w);
02109 gauss_points_tr (gp1.a, gp2.a, w.a, intordsm[ii][jj]);
02110 for (k = 0; k < intordsm[ii][jj]; k++)
02111 {
02112 xi=gp1[k];
02113 eta=gp2[k];
02114
02115 ipval = approx_nat (xi,eta,anv);
02116 ncompstr = Mm->ip[ipp].ncompstr;
02117 ncompeqother = Mm->ip[ipp].ncompeqother;
02118 if ((ictn[0] & inistrain) && (j < ncompstr))
02119 {
02120 Mm->ip[ipp].strain[idstra] += ipval;
02121 ipp++;
02122 continue;
02123 }
02124 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02125 {
02126 Mm->ip[ipp].stress[idstre] += ipval;
02127 ipp++;
02128 continue;
02129 }
02130 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
02131 {
02132 Mm->ip[ipp].eqother[idoth] += ipval;
02133 ipp++;
02134 continue;
02135 }
02136 if ((ictn[0] & inicond) && (j < nv))
02137 {
02138 if (Mm->ic[ipp] == NULL)
02139 {
02140 Mm->ic[ipp] = new double[nv-j];
02141 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
02142 }
02143 Mm->ic[ipp][idic] += ipval;
02144 ipp++;
02145 continue;
02146 }
02147 ipp++;
02148 }
02149 destrv(gp1); destrv (gp2); destrv (w);
02150 }
02151 }
02152 ipp=Mt->elements[eid].ipp[ri][ci];
02153 ncompstr = Mm->ip[ipp].ncompstr;
02154 ncompeqother = Mm->ip[ipp].ncompeqother;
02155 if ((ictn[0] & inistrain) && (j < ncompstr))
02156 {
02157 nstra++;
02158 idstra++;
02159 continue;
02160 }
02161 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02162 {
02163 nstre++;
02164 idstre++;
02165 continue;
02166 }
02167 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
02168 {
02169 idoth++;
02170 continue;
02171 }
02172 if ((ictn[0] & inicond) && (j < nv))
02173 {
02174 idic++;
02175 continue;
02176 }
02177 }
02178 }
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188 void planeelemrotlt::ipvolume (long eid,long ri,long ci)
02189 {
02190 long i,ii,jj,ipp;
02191 double jac,thick;
02192 ivector nodes(nne);
02193 vector x(nne),y(nne),t(nne),gp1,gp2,w;
02194
02195
02196 Mt->give_elemnodes (eid,nodes);
02197
02198 Mc->give_thickness (eid,nodes,t);
02199 Mt->give_node_coord2d (x,y,eid);
02200
02201 for (ii=0;ii<nb;ii++){
02202 for (jj=0;jj<nb;jj++){
02203 if (intordsm[ii][jj]==0) continue;
02204
02205 allocv (nip[ii][jj],gp1);
02206 allocv (nip[ii][jj],gp2);
02207 allocv (nip[ii][jj],w);
02208
02209
02210 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
02211 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
02212
02213 for (i=0;i<intordsm[ii][jj];i++){
02214 thick = approx_nat (gp1[i],gp2[i],t);
02215 jac_2d (jac,x,y,gp1[i],gp2[i]);
02216 jac*=w[i]*thick;
02217
02218 Mm->storeipvol (ipp,jac);
02219
02220 ipp++;
02221 }
02222
02223 destrv (w); destrv (gp1); destrv (gp2);
02224 }
02225 }
02226
02227 allocv (1,w);
02228 allocv (1,gp1);
02229 allocv (1,gp2);
02230
02231 gauss_points_tr (gp1.a,gp2.a,w.a,1);
02232
02233 thick = approx_nat (gp1[0],gp2[0],t);
02234 jac_2d (jac,x,y,gp1[0],gp2[0]);
02235 jac*=w[0]*thick;
02236
02237 Mm->storeipvol (ipp,jac);
02238
02239 destrv (gp1); destrv (gp2); destrv (w);
02240 }