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