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