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