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