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