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