00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include "linhex.h"
00004 #include "global.h"
00005 #include "globmat.h"
00006 #include "genfile.h"
00007 #include "intpoints.h"
00008 #include "node.h"
00009 #include "element.h"
00010 #include "loadcase.h"
00011 #include "matrix.h"
00012
00013
00014
00015
00016 linhex::linhex ()
00017 {
00018 long i,j;
00019
00020
00021 nne=8;
00022
00023 ndofe=24;
00024
00025 tncomp=6;
00026
00027 napfun=3;
00028
00029 intordmm=2;
00030
00031 ned=12;
00032
00033 nned=2;
00034
00035 intordb=2;
00036
00037 nsurf=6;
00038
00039 nnsurf=4;
00040
00041 ssst=spacestress;
00042
00043
00044 nb=1;
00045
00046
00047 ncomp = new long [nb];
00048 ncomp[0]=6;
00049
00050
00051 cncomp = new long [nb];
00052 cncomp[0]=0;
00053
00054
00055
00056 nip = new long* [nb];
00057 intordsm = new long* [nb];
00058 for (i=0;i<nb;i++){
00059 nip[i] = new long [nb];
00060 intordsm[i] = new long [nb];
00061 }
00062
00063 nip[0][0]=8;
00064
00065
00066 tnip=0;
00067 for (i=0;i<nb;i++){
00068 for (j=0;j<nb;j++){
00069 tnip+=nip[i][j];
00070 }
00071 }
00072
00073 intordsm[0][0]=2;
00074 }
00075
00076 linhex::~linhex ()
00077 {
00078 long i;
00079
00080 for (i=0;i<nb;i++){
00081 delete [] nip[i];
00082 delete [] intordsm[i];
00083 }
00084 delete [] nip;
00085 delete [] intordsm;
00086
00087 delete [] cncomp;
00088 delete [] ncomp;
00089 }
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 double linhex::approx (double xi,double eta,double zeta,vector &nodval)
00101 {
00102 double f;
00103 vector bf(nne);
00104
00105 bf_lin_hex_3d (bf.a,xi,eta,zeta);
00106
00107 scprd (bf,nodval,f);
00108
00109 return f;
00110 }
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 void linhex::bf_matrix (matrix &n,double xi,double eta,double zeta)
00121 {
00122 long i,j,k,l;
00123 vector bf(nne);
00124
00125 fillm (0.0,n);
00126
00127 bf_lin_hex_3d (bf.a,xi,eta,zeta);
00128
00129 j=0; k=1; l=2;
00130 for (i=0;i<nne;i++){
00131 n[0][j]=bf[i]; j+=3;
00132 n[1][k]=bf[i]; k+=3;
00133 n[2][l]=bf[i]; l+=3;
00134 }
00135 }
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147 void linhex::geom_matrix (matrix &gm,vector &x,vector &y,vector &z,
00148 double xi,double eta,double zeta,double &jac)
00149 {
00150 long i,j,k,l;
00151 vector dx(nne),dy(nne),dz(nne);
00152
00153 dx_bf_lin_hex_3d (dx.a,eta,zeta);
00154 dy_bf_lin_hex_3d (dy.a,xi,zeta);
00155 dz_bf_lin_hex_3d (dz.a,xi,eta);
00156
00157 derivatives_3d (dx,dy,dz,jac,x,y,z,xi,eta,zeta);
00158
00159 fillm (0.0,gm);
00160
00161 j=0; k=1; l=2;
00162 for (i=0;i<nne;i++){
00163 gm[0][j]=dx[i];
00164 gm[1][k]=dy[i];
00165 gm[2][l]=dz[i];
00166
00167 gm[3][k]=dz[i];
00168 gm[3][l]=dy[i];
00169
00170 gm[4][j]=dz[i];
00171 gm[4][l]=dx[i];
00172
00173 gm[5][j]=dy[i];
00174 gm[5][k]=dx[i];
00175
00176 j+=3; k+=3; l+=3;
00177 }
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 void linhex::bvectors (vector &x,vector &y,vector &z,double xi,double eta,double zeta,double &jac,
00192 vector &b11,vector &b12,vector &b13,
00193 vector &b21,vector &b22,vector &b23,
00194 vector &b31,vector &b32,vector &b33)
00195 {
00196 vector dx(nne),dy(nne),dz(nne);
00197
00198 dx_bf_lin_hex_3d (dx.a,eta,zeta);
00199 dy_bf_lin_hex_3d (dy.a,xi,zeta);
00200 dz_bf_lin_hex_3d (dz.a,xi,eta);
00201
00202 derivatives_3d (dx,dy,dz,jac,x,y,z,xi,eta,zeta);
00203
00204 fillv (0.0,b11);
00205 fillv (0.0,b12);
00206 fillv (0.0,b13);
00207 fillv (0.0,b21);
00208 fillv (0.0,b22);
00209 fillv (0.0,b23);
00210 fillv (0.0,b31);
00211 fillv (0.0,b32);
00212 fillv (0.0,b33);
00213
00214
00215 b11[0]=dx[0]; b11[3]=dx[1]; b11[6]=dx[2]; b11[9]=dx[3]; b11[12]=dx[4]; b11[15]=dx[5]; b11[18]=dx[6]; b11[21]=dx[7];
00216
00217 b12[0]=dy[0]; b12[3]=dy[1]; b12[6]=dy[2]; b12[9]=dy[3]; b12[12]=dy[4]; b12[15]=dy[5]; b12[18]=dy[6]; b12[21]=dy[7];
00218
00219 b13[0]=dz[0]; b13[3]=dz[1]; b13[6]=dz[2]; b13[9]=dz[3]; b13[12]=dz[4]; b13[15]=dz[5]; b13[18]=dz[6]; b13[21]=dz[7];
00220
00221
00222 b21[1]=dx[0]; b21[4]=dx[1]; b21[7]=dx[2]; b21[10]=dx[3]; b21[13]=dx[4]; b21[16]=dx[5]; b21[19]=dx[6]; b21[22]=dx[7];
00223
00224 b22[1]=dy[0]; b22[4]=dy[1]; b22[7]=dy[2]; b22[10]=dy[3]; b22[13]=dy[4]; b22[16]=dy[5]; b22[19]=dy[6]; b22[22]=dy[7];
00225
00226 b23[1]=dz[0]; b23[4]=dz[1]; b23[7]=dz[2]; b23[10]=dz[3]; b23[13]=dz[4]; b23[16]=dz[5]; b23[19]=dz[6]; b23[22]=dz[7];
00227
00228
00229 b31[2]=dx[0]; b31[5]=dx[1]; b31[8]=dx[2]; b31[11]=dx[3]; b31[14]=dx[4]; b31[17]=dx[5]; b31[20]=dx[6]; b31[23]=dx[7];
00230
00231 b32[2]=dy[0]; b32[5]=dy[1]; b32[8]=dy[2]; b32[11]=dy[3]; b32[14]=dy[4]; b32[17]=dy[5]; b32[20]=dy[6]; b32[23]=dy[7];
00232
00233 b33[2]=dz[0]; b33[5]=dz[1]; b33[8]=dz[2]; b33[11]=dz[3]; b33[14]=dz[4]; b33[17]=dz[5]; b33[20]=dz[6]; b33[23]=dz[7];
00234
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 void linhex::gngeom_matrix (matrix &gm,vector &r,vector &x,vector &y,vector &z,double xi,double eta,double zeta,double &jac)
00249 {
00250 long i;
00251 double b11r,b12r,b13r,b21r,b22r,b23r,b31r,b32r,b33r;
00252 vector b11(ndofe),b12(ndofe),b13(ndofe),b21(ndofe),b22(ndofe),b23(ndofe),b31(ndofe),b32(ndofe),b33(ndofe),av(ndofe);
00253
00254 fillm (0.0,gm);
00255
00256 bvectors (x,y,z,xi,eta,zeta,jac,b11,b12,b13,b21,b22,b23,b31,b32,b33);
00257
00258 scprd (b11,r,b11r);
00259 scprd (b12,r,b12r);
00260 scprd (b13,r,b13r);
00261 scprd (b21,r,b21r);
00262 scprd (b22,r,b22r);
00263 scprd (b23,r,b23r);
00264 scprd (b31,r,b31r);
00265 scprd (b32,r,b32r);
00266 scprd (b33,r,b33r);
00267
00268
00269
00270
00271
00272
00273
00274 for (i=0;i<ndofe;i++){
00275 gm[0][i]+=b11[i];
00276 }
00277
00278
00279 cmulv(b11r,b11,av);
00280 for (i=0;i<ndofe;i++){
00281 gm[0][i]+=av[i];
00282 }
00283
00284
00285 cmulv(b21r,b21,av);
00286 for (i=0;i<ndofe;i++){
00287 gm[0][i]+=av[i];
00288 }
00289
00290
00291 cmulv(b31r,b31,av);
00292 for (i=0;i<ndofe;i++){
00293 gm[0][i]+=av[i];
00294 }
00295
00296
00297
00298
00299
00300
00301
00302 for (i=0;i<ndofe;i++){
00303 gm[1][i]+=b22[i];
00304 }
00305
00306
00307 cmulv(b12r,b12,av);
00308 for (i=0;i<ndofe;i++){
00309 gm[1][i]+=av[i];
00310 }
00311
00312
00313 cmulv(b22r,b22,av);
00314 for (i=0;i<ndofe;i++){
00315 gm[1][i]+=av[i];
00316 }
00317
00318
00319 cmulv(b32r,b32,av);
00320 for (i=0;i<ndofe;i++){
00321 gm[1][i]+=av[i];
00322 }
00323
00324
00325
00326
00327
00328
00329 for (i=0;i<ndofe;i++){
00330 gm[2][i]+=b22[i];
00331 }
00332
00333
00334 cmulv(b13r,b13,av);
00335 for (i=0;i<ndofe;i++){
00336 gm[2][i]+=av[i];
00337 }
00338
00339
00340 cmulv(b23r,b23,av);
00341 for (i=0;i<ndofe;i++){
00342 gm[2][i]+=av[i];
00343 }
00344
00345
00346 cmulv(b33r,b33,av);
00347 for (i=0;i<ndofe;i++){
00348 gm[2][i]+=av[i];
00349 }
00350
00351
00352
00353
00354
00355
00356
00357 for (i=0;i<ndofe;i++){
00358 gm[3][i]+=b23[i]+b32[i];
00359 }
00360
00361
00362 cmulv(b13r,b12,av);
00363 for (i=0;i<ndofe;i++){
00364 gm[3][i]+=av[i];
00365 }
00366
00367
00368 cmulv(b12r,b13,av);
00369 for (i=0;i<ndofe;i++){
00370 gm[3][i]+=av[i];
00371 }
00372
00373
00374 cmulv(b23r,b22,av);
00375 for (i=0;i<ndofe;i++){
00376 gm[3][i]+=av[i];
00377 }
00378
00379
00380 cmulv(b22r,b23,av);
00381 for (i=0;i<ndofe;i++){
00382 gm[3][i]+=av[i];
00383 }
00384
00385
00386 cmulv(b33r,b32,av);
00387 for (i=0;i<ndofe;i++){
00388 gm[3][i]+=av[i];
00389 }
00390
00391
00392 cmulv(b32r,b33,av);
00393 for (i=0;i<ndofe;i++){
00394 gm[3][i]+=av[i];
00395 }
00396
00397
00398
00399
00400
00401
00402 for (i=0;i<ndofe;i++){
00403 gm[4][i]+=b31[i]+b13[i];
00404 }
00405
00406
00407 cmulv(b11r,b13,av);
00408 for (i=0;i<ndofe;i++){
00409 gm[4][i]+=av[i];
00410 }
00411
00412
00413 cmulv(b13r,b11,av);
00414 for (i=0;i<ndofe;i++){
00415 gm[4][i]+=av[i];
00416 }
00417
00418
00419 cmulv(b21r,b23,av);
00420 for (i=0;i<ndofe;i++){
00421 gm[4][i]+=av[i];
00422 }
00423
00424
00425 cmulv(b23r,b21,av);
00426 for (i=0;i<ndofe;i++){
00427 gm[4][i]+=av[i];
00428 }
00429
00430
00431 cmulv(b31r,b33,av);
00432 for (i=0;i<ndofe;i++){
00433 gm[4][i]+=av[i];
00434 }
00435
00436
00437 cmulv(b33r,b31,av);
00438 for (i=0;i<ndofe;i++){
00439 gm[4][i]+=av[i];
00440 }
00441
00442
00443
00444
00445
00446
00447
00448 for (i=0;i<ndofe;i++){
00449 gm[5][i]+=b12[i]+b21[i];
00450 }
00451
00452
00453 cmulv(b12r,b11,av);
00454 for (i=0;i<ndofe;i++){
00455 gm[5][i]+=av[i];
00456 }
00457
00458
00459 cmulv(b11r,b12,av);
00460 for (i=0;i<ndofe;i++){
00461 gm[5][i]+=av[i];
00462 }
00463
00464
00465 cmulv(b22r,b21,av);
00466 for (i=0;i<ndofe;i++){
00467 gm[5][i]+=av[i];
00468 }
00469
00470
00471 cmulv(b21r,b22,av);
00472 for (i=0;i<ndofe;i++){
00473 gm[5][i]+=av[i];
00474 }
00475
00476
00477 cmulv(b32r,b31,av);
00478 for (i=0;i<ndofe;i++){
00479 gm[5][i]+=av[i];
00480 }
00481
00482
00483 cmulv(b31r,b32,av);
00484 for (i=0;i<ndofe;i++){
00485 gm[5][i]+=av[i];
00486 }
00487
00488
00489
00490 }
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503 void linhex::gnl_grmatrix (matrix &grm,vector &x,vector &y,vector &z,double xi,double eta,double zeta,double &jac)
00504 {
00505 long i;
00506 vector b11(ndofe),b12(ndofe),b13(ndofe),b21(ndofe),b22(ndofe),b23(ndofe),b31(ndofe),b32(ndofe),b33(ndofe);
00507
00508 bvectors (x,y,z,xi,eta,zeta,jac,b11,b12,b13,b21,b22,b23,b31,b32,b33);
00509
00510 for (i=0;i<ndofe;i++){
00511 grm[0][i]=b11[i];
00512 grm[1][i]=b12[i];
00513 grm[2][i]=b13[i];
00514 grm[3][i]=b21[i];
00515 grm[4][i]=b22[i];
00516 grm[5][i]=b23[i];
00517 grm[6][i]=b31[i];
00518 grm[7][i]=b32[i];
00519 grm[8][i]=b33[i];
00520 }
00521 }
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535 void linhex::transf_matrix (ivector &nodes,matrix &tmat)
00536 {
00537 long i,n,m;
00538
00539 fillm (0.0,tmat);
00540
00541 n=nodes.n;
00542 m=tmat.m;
00543 for (i=0;i<m;i++){
00544 tmat[i][i]=1.0;
00545 }
00546
00547 for (i=0;i<n;i++){
00548 if (Mt->nodes[nodes[i]].transf>0){
00549 tmat[i*3+0][i*3]=Mt->nodes[nodes[i]].e1[0];
00550 tmat[i*3+1][i*3]=Mt->nodes[nodes[i]].e1[1];
00551 tmat[i*3+2][i*3]=Mt->nodes[nodes[i]].e1[2];
00552
00553 tmat[i*3+0][i*3+1]=Mt->nodes[nodes[i]].e2[0];
00554 tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e2[1];
00555 tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e2[2];
00556
00557 tmat[i*3+0][i*3+2]=Mt->nodes[nodes[i]].e3[0];
00558 tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e3[1];
00559 tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e3[2];
00560 }
00561 }
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575 void linhex::gl_stiffness_matrix (long eid,long ri,long ci,matrix &sm)
00576 {
00577 long i,j,k,ii,jj,ipp;
00578 double xi,eta,zeta,jac;
00579 vector x(nne),y(nne),z(nne),w,gp;
00580 matrix gm,d(tncomp,tncomp);
00581
00582 Mt->give_node_coord3d (x,y,z,eid);
00583
00584 fillm (0.0,sm);
00585
00586 for (ii=0;ii<nb;ii++){
00587 allocm (ncomp[ii],ndofe,gm);
00588 for (jj=0;jj<nb;jj++){
00589 if (intordsm[ii][jj]==0) continue;
00590
00591 allocv (intordsm[ii][jj],w);
00592 allocv (intordsm[ii][jj],gp);
00593
00594 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00595
00596 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00597
00598 for (i=0;i<intordsm[ii][jj];i++){
00599 xi=gp[i];
00600 for (j=0;j<intordsm[ii][jj];j++){
00601 eta=gp[j];
00602 for (k=0;k<intordsm[ii][jj];k++){
00603 zeta=gp[k];
00604
00605
00606 geom_matrix (gm,x,y,z,xi,eta,zeta,jac);
00607
00608
00609 Mm->matstiff (d,ipp); ipp++;
00610
00611 jac*=w[i]*w[j]*w[k];
00612
00613 if (jac<0.0){
00614
00615 }
00616
00617
00618 bdbjac (sm,gm,d,gm,jac);
00619
00620 }
00621 }
00622 }
00623 destrv (gp); destrv (w);
00624 }
00625 destrm (gm);
00626 }
00627 }
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641 void linhex::gnl_stiffness_matrix (long lcid,long eid,long ri,long ci,matrix &sm)
00642 {
00643 long i,j,k,ipp;
00644 double xi,eta,zeta,jac,jac2;
00645 vector w,gp,sig(tncomp),r(ndofe),x(nne),y(nne),z(nne);
00646 matrix gm(tncomp,ndofe),grm(9,ndofe),d(tncomp,tncomp),s(9,9);
00647
00648
00649 Mt->give_node_coord3d (x,y,z,eid);
00650
00651 eldispl (lcid,eid,r.a);
00652
00653
00654 fillm (0.0,sm);
00655
00656
00657 allocv (intordsm[0][0],w);
00658
00659 allocv (intordsm[0][0],gp);
00660
00661
00662 gauss_points (gp.a,w.a,intordsm[0][0]);
00663
00664
00665 ipp=Mt->elements[eid].ipp[ri][ci];
00666
00667
00668 for (i=0;i<intordsm[0][0];i++){
00669 xi=gp[i];
00670 for (j=0;j<intordsm[0][0];j++){
00671 eta=gp[j];
00672 for (k=0;k<intordsm[0][0];k++){
00673 zeta=gp[k];
00674
00675
00676
00677
00678
00679
00680 gngeom_matrix (gm,r,x,y,z,xi,eta,zeta,jac);
00681
00682
00683 Mm->matstiff (d,ipp);
00684
00685 jac*=w[i]*w[j]*w[k];
00686
00687
00688 bdbjac (sm,gm,d,gm,jac);
00689
00690
00691
00692
00693
00694
00695
00696 gnl_grmatrix (grm,x,y,z,xi,eta,zeta,jac2);
00697
00698
00699 Mm->givestress (lcid,ipp,sig);
00700
00701 s[0][0]=sig[0]; s[0][1]=sig[5]; s[0][2]=sig[4];
00702 s[1][0]=sig[5]; s[1][1]=sig[1]; s[1][2]=sig[3];
00703 s[2][0]=sig[4]; s[2][1]=sig[3]; s[2][2]=sig[2];
00704
00705 s[3][3]=sig[0]; s[3][4]=sig[5]; s[3][5]=sig[4];
00706 s[4][3]=sig[5]; s[4][4]=sig[1]; s[4][5]=sig[3];
00707 s[5][3]=sig[4]; s[5][4]=sig[3]; s[5][5]=sig[2];
00708
00709 s[6][6]=sig[0]; s[6][7]=sig[5]; s[6][8]=sig[4];
00710 s[7][6]=sig[5]; s[7][7]=sig[1]; s[7][8]=sig[3];
00711 s[8][6]=sig[4]; s[8][7]=sig[3]; s[8][8]=sig[2];
00712
00713
00714
00715
00716 bdbjac (sm,grm,s,grm,jac);
00717
00718 ipp++;
00719 }
00720 }
00721 }
00722
00723 destrv (gp); destrv (w);
00724
00725 }
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737 void linhex::res_stiffness_matrix (long lcid,long eid,matrix &sm)
00738 {
00739 gl_stiffness_matrix (eid,0,0,sm);
00740
00741 ivector nodes (nne);
00742 Mt->give_elemnodes (eid,nodes);
00743 long transf = Mt->locsystems (nodes);
00744 if (transf>0){
00745 matrix tmat (ndofe,ndofe);
00746 transf_matrix (nodes,tmat);
00747 glmatrixtransf (sm,tmat);
00748 }
00749
00750 }
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760 void linhex::mass_matrix (long eid,matrix &mm)
00761 {
00762 long i,j,k;
00763 double jac,xi,eta,zeta,rho;
00764 ivector nodes (nne);
00765 vector x(nne),y(nne),z(nne),w(intordmm),gp(intordmm),dens(nne);
00766 matrix n(napfun,ndofe);
00767
00768 Mt->give_elemnodes (eid,nodes);
00769 Mc->give_density (eid,nodes,dens);
00770
00771 Mt->give_node_coord3d (x,y,z,eid);
00772 gauss_points (gp.a,w.a,intordmm);
00773
00774 fillm (0.0,mm);
00775
00776 for (i=0;i<intordmm;i++){
00777 xi=gp[i];
00778 for (j=0;j<intordmm;j++){
00779 eta=gp[j];
00780 for (k=0;k<intordmm;k++){
00781 zeta=gp[k];
00782
00783 jac_3d (jac,x,y,z,xi,eta,zeta);
00784
00785 jac=fabs(jac);
00786
00787 bf_matrix (n,xi,eta,zeta);
00788
00789 rho = approx (xi,eta,zeta,dens);
00790
00791 jac*=w[i]*w[j]*w[k]*rho;
00792
00793 nnj (mm.a,n.a,jac,n.m,n.n);
00794 }
00795 }
00796 }
00797
00798 }
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808 void linhex::res_mass_matrix (long eid,matrix &mm)
00809 {
00810 long transf;
00811
00812 mass_matrix (eid,mm);
00813
00814
00815
00816 ivector nodes (nne);
00817 Mt->give_elemnodes (eid,nodes);
00818 transf = Mt->locsystems (nodes);
00819 if (transf>0){
00820 matrix tmat (ndofe,ndofe);
00821 transf_matrix (nodes,tmat);
00822 glmatrixtransf (mm,tmat);
00823 }
00824 }
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834 void linhex::load_matrix (long eid,matrix &lm)
00835 {
00836 long i,j,k;
00837 double jac,xi,eta,zeta,w1,w2,w3;
00838 ivector nodes (nne);
00839 vector x(nne),y(nne),z(nne),w(intordmm),gp(intordmm);
00840 matrix n(napfun,ndofe);
00841
00842 Mt->give_elemnodes (eid,nodes);
00843 Mt->give_node_coord3d (x,y,z,eid);
00844 gauss_points (gp.a,w.a,intordmm);
00845
00846 fillm (0.0,lm);
00847
00848 for (i=0;i<intordmm;i++){
00849 xi=gp[i]; w1=w[i];
00850 for (j=0;j<intordmm;j++){
00851 eta=gp[j]; w2=w[j];
00852 for (k=0;k<intordmm;k++){
00853 zeta=gp[k]; w3=w[k];
00854
00855 jac_3d (jac,x,y,z,xi,eta,zeta);
00856
00857 bf_matrix (n,xi,eta,zeta);
00858
00859 jac*=w1*w2*w3;
00860
00861 nnj (lm.a,n.a,jac,n.m,n.n);
00862 }
00863 }
00864 }
00865
00866 }
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876 void linhex::res_load_matrix (long eid,matrix &lm)
00877 {
00878 long transf;
00879 ivector nodes(nne);
00880
00881 load_matrix (eid,lm);
00882
00883
00884
00885 Mt->give_elemnodes (eid,nodes);
00886 transf = Mt->locsystems (nodes);
00887 if (transf>0){
00888 matrix tmat (ndofe,ndofe);
00889 transf_matrix (nodes,tmat);
00890 glmatrixtransf (lm,tmat);
00891 }
00892 }
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905 void linhex::res_ip_strains (long lcid,long eid)
00906 {
00907 ivector nodes(nne);
00908 vector x(nne),y(nne),z(nne),r(ndofe),aux;
00909 matrix tmat;
00910
00911 Mt->give_elemnodes (eid,nodes);
00912 Mt->give_node_coord3d (x,y,z,eid);
00913 eldispl (lcid,eid,r.a);
00914
00915
00916
00917 long transf = Mt->locsystems (nodes);
00918 if (transf>0){
00919 allocv (ndofe,aux);
00920 allocm (ndofe,ndofe,tmat);
00921 transf_matrix (nodes,tmat);
00922
00923 lgvectortransf (aux,r,tmat);
00924 copyv (aux,r);
00925 destrv (aux);
00926 destrm (tmat);
00927 }
00928
00929 gl_ip_strains (lcid,eid,0,0,x,y,z,r);
00930
00931
00932 }
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947 void linhex::gl_ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &z,vector &r)
00948 {
00949 long i,j,k,ii,ipp;
00950 double xi,eta,zeta,jac;
00951 vector gp,w,eps;
00952 matrix gm;
00953
00954 for (ii=0;ii<nb;ii++){
00955
00956 allocv (intordsm[ii][ii],gp);
00957 allocv (intordsm[ii][ii],w);
00958 allocv (ncomp[ii],eps);
00959 allocm (ncomp[ii],ndofe,gm);
00960
00961 gauss_points (gp.a,w.a,intordsm[ii][ii]);
00962
00963 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00964 for (i=0;i<intordsm[ii][ii];i++){
00965 xi=gp[i];
00966 for (j=0;j<intordsm[ii][ii];j++){
00967 eta=gp[j];
00968 for (k=0;k<intordsm[ii][ii];k++){
00969 zeta=gp[k];
00970
00971 geom_matrix (gm,x,y,z,xi,eta,zeta,jac);
00972
00973 mxv (gm,r,eps);
00974
00975 Mm->storestrain (lcid,ipp,cncomp[ii],eps);
00976
00977 ipp++;
00978 }
00979 }
00980 }
00981
00982 destrm (gm); destrv (w); destrv (gp);
00983 destrv (eps);
00984 }
00985 }
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000 void linhex::gnl_ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &z,vector &r)
01001 {
01002 long i,j,k,ipp;
01003 double xi,eta,zeta,jac,b11r,b12r,b13r,b21r,b22r,b23r,b31r,b32r,b33r;
01004 vector gp,w,eps;
01005 vector b11(ndofe),b12(ndofe),b13(ndofe),b21(ndofe),b22(ndofe),b23(ndofe),b31(ndofe),b32(ndofe),b33(ndofe);
01006
01007 allocv (intordsm[0][0],gp);
01008 allocv (intordsm[0][0],w);
01009 allocv (tncomp,eps);
01010
01011 gauss_points (gp.a,w.a,intordsm[0][0]);
01012
01013 ipp=Mt->elements[eid].ipp[ri][ci];
01014 for (i=0;i<intordsm[0][0];i++){
01015 xi=gp[i];
01016 for (j=0;j<intordsm[0][0];j++){
01017 eta=gp[j];
01018 for (k=0;k<intordsm[0][0];k++){
01019 zeta=gp[k];
01020
01021 bvectors (x,y,z,xi,eta,zeta,jac,b11,b12,b13,b21,b22,b23,b31,b32,b33);
01022
01023 scprd (b11,r,b11r);
01024 scprd (b12,r,b12r);
01025 scprd (b13,r,b13r);
01026 scprd (b21,r,b21r);
01027 scprd (b22,r,b22r);
01028 scprd (b23,r,b23r);
01029 scprd (b31,r,b31r);
01030 scprd (b32,r,b32r);
01031 scprd (b33,r,b33r);
01032
01033 eps[0] = b11r + 0.5*b11r*b11r + 0.5*b21r*b21r + 0.5*b31r*b31r;
01034 eps[1] = b22r + 0.5*b12r*b12r + 0.5*b22r*b22r + 0.5*b32r*b32r;
01035 eps[2] = b33r + 0.5*b13r*b13r + 0.5*b23r*b23r + 0.5*b33r*b33r;
01036
01037 eps[3] = b23r+b32r + b12r*b13r + b22r*b23r + b32r*b33r;
01038 eps[4] = b31r+b13r + b13r*b11r + b23r*b21r + b33r*b31r;
01039 eps[5] = b12r+b21r + b11r*b12r + b21r*b22r + b31r*b32r;
01040
01041
01042 Mm->storestrain (lcid,ipp,eps);
01043 ipp++;
01044 }
01045 }
01046 }
01047
01048 destrv (eps); destrv (w); destrv (gp);
01049
01050 }
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062 void linhex::nod_strains_ip (long lcid,long eid,long ri,long ci)
01063 {
01064 long i,j,ipp;
01065 ivector ipnum(nne),nod(nne);
01066 vector eps(tncomp);
01067
01068
01069
01070 ipp=Mt->elements[eid].ipp[ri][ci];
01071 nodip_linhex (ipp,intordsm[0][0],ipnum);
01072
01073
01074 Mt->give_elemnodes (eid,nod);
01075
01076 for (i=0;i<nne;i++){
01077
01078 Mm->givestrain (lcid,ipnum[i],eps);
01079
01080
01081 j=nod[i];
01082 Mt->nodes[j].storestrain (lcid,0,eps);
01083 }
01084 }
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097 void linhex::nod_strains_comp (long lcid,long eid,double **stra)
01098 {
01099 long i,j;
01100 double jac;
01101 ivector nodes(nne);
01102 vector x(nne),y(nne),z(nne),nxi(nne),neta(nne),nzeta(nne),r(ndofe),eps(tncomp),aux;
01103 matrix tmat,gm(tncomp,ndofe);
01104
01105
01106 Mt->give_node_coord3d (x,y,z,eid);
01107
01108 Mt->give_elemnodes (eid,nodes);
01109
01110 eldispl (lcid,eid,r.a);
01111
01112
01113 long transf = Mt->locsystems (nodes);
01114 if (transf>0){
01115 allocv (ndofe,aux);
01116 allocm (ndofe,ndofe,tmat);
01117 transf_matrix (nodes,tmat);
01118
01119 lgvectortransf (aux,r,tmat);
01120 copyv (aux,r);
01121 destrv (aux);
01122 destrm (tmat);
01123 }
01124
01125
01126
01127 nodcoord_linhex (nxi,neta,nzeta);
01128
01129
01130 for (i=0;i<nne;i++){
01131
01132 geom_matrix (gm,x,y,z,nxi[i],neta[i],nzeta[i],jac);
01133
01134 mxv (gm,r,eps);
01135
01136 for (j=0;j<eps.n;j++){
01137 stra[i][j]=eps[j];
01138 }
01139 }
01140
01141 }
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157 void linhex::appstrain (long lcid,long eid,double xi,double eta,double zeta,long fi,long ncomp,vector &eps)
01158 {
01159 long i,j,k;
01160 ivector nod(nne);
01161 vector nodval(nne);
01162
01163 if (ncomp != eps.n){
01164 fprintf (stderr,"\n\n wrong interval of indices in function strain (%s, line %d).\n",__FILE__,__LINE__);
01165 abort ();
01166 }
01167
01168 Mt->give_elemnodes (eid,nod);
01169 k=0;
01170 for (i=fi;i<fi+ncomp;i++){
01171 for (j=0;j<nne;j++){
01172 nodval[j]=Mt->nodes[nod[j]].strain[lcid*tncomp+i];
01173 }
01174 eps[k]=approx (xi,eta,zeta,nodval);
01175 k++;
01176 }
01177 }
01178
01179
01180
01181 void linhex::strains (long lcid,long eid,long ri,long ci)
01182 {
01183 long i,naep,ncp,sid;
01184 double **stra;
01185 vector coord,eps;
01186
01187 if (Mp->strainaver==0){
01188 stra = new double* [nne];
01189 for (i=0;i<nne;i++){
01190 stra[i] = new double [tncomp];
01191 }
01192
01193 }
01194
01195 switch (Mm->stra.tape[eid]){
01196 case nowhere:{
01197 break;
01198 }
01199 case intpts:{
01200
01201 break;
01202 }
01203 case enodes:{
01204 break;
01205 }
01206 case userdefined:{
01207
01208 naep = Mm->stra.give_naep (eid);
01209 ncp = Mm->stra.give_ncomp (eid);
01210 sid = Mm->stra.give_sid (eid);
01211 allocv (ncp,eps);
01212 allocv (3,coord);
01213 for (i=0;i<naep;i++){
01214 Mm->stra.give_aepcoord (sid,i,coord);
01215
01216 if (Mp->strainaver==0)
01217
01218 if (Mp->strainaver==1)
01219 appstrain (lcid,eid,coord[0],coord[1],coord[2],0,ncp,eps);
01220
01221 Mm->stra.storevalues(lcid,eid,i,eps);
01222 }
01223 destrv (eps);
01224 destrv (coord);
01225 break;
01226 }
01227 default:{
01228 fprintf (stderr,"\n\n unknown strain point is required in function planeelemlq::strains (%s, line %d).\n",__FILE__,__LINE__);
01229 }
01230 }
01231
01232 if (Mp->strainaver==0){
01233 for (i=0;i<nne;i++){
01234 delete [] stra[i];
01235 }
01236 delete [] stra;
01237 }
01238 }
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248 void linhex::res_ip_stresses (long lcid,long eid)
01249 {
01250 ip_stresses (lcid,eid,0,0);
01251 }
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264 void linhex::ip_stresses (long lcid,long eid,long ri,long ci)
01265 {
01266 long i,j,k,ipp;
01267
01268 ipp=Mt->elements[eid].ipp[ri][ci];
01269
01270 for (i=0;i<intordsm[0][0];i++){
01271 for (j=0;j<intordsm[0][0];j++){
01272 for (k=0;k<intordsm[0][0];k++){
01273
01274
01275 if (Mp->strcomp==1)
01276 Mm->computenlstresses (ipp);
01277
01278 ipp++;
01279 }
01280 }
01281 }
01282 }
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295 void linhex::ip_elast_stresses (long lcid,long eid,long ri,long ci)
01296 {
01297 long i,j,k,ipp;
01298 vector gp,w,eps(tncomp),sig(tncomp);
01299 matrix d(tncomp,tncomp);
01300
01301 allocv (intordsm[0][0],gp);
01302 allocv (intordsm[0][0],w);
01303
01304 gauss_points (gp.a,w.a,intordsm[0][0]);
01305 ipp=Mt->elements[eid].ipp[ri][ci];
01306
01307 for (i=0;i<intordsm[0][0];i++){
01308 for (j=0;j<intordsm[0][0];j++){
01309 for (k=0;k<intordsm[0][0];k++){
01310
01311 Mm->matstiff (d,ipp);
01312
01313 Mm->givestrain (lcid,ipp,eps);
01314
01315 mxv (d,eps,sig);
01316
01317 Mm->storestress (lcid,ipp,sig);
01318
01319 ipp++;
01320 }
01321 }
01322 }
01323
01324 destrv (w); destrv (gp);
01325 }
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336 void linhex::nod_stresses_ip (long lcid,long eid,long ri,long ci)
01337 {
01338 long i,j,ipp;
01339 ivector ipnum(nne),nod(nne);
01340 vector sig(tncomp);
01341
01342
01343
01344 ipp=Mt->elements[eid].ipp[ri][ci];
01345 nodip_linhex (ipp,intordsm[0][0],ipnum);
01346
01347
01348 Mt->give_elemnodes (eid,nod);
01349
01350 for (i=0;i<nne;i++){
01351
01352 Mm->givestress (lcid,ipnum[i],sig);
01353
01354
01355 j=nod[i];
01356 Mt->nodes[j].storestress (lcid,0,sig);
01357 }
01358 }
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372 void linhex::appstress (long lcid,long eid,double xi,double eta,double zeta,long fi,long ncomp,vector &sig)
01373 {
01374 long i,j,k;
01375 ivector nodes(nne);
01376 vector nodval(nne);
01377
01378 if (ncomp != sig.n){
01379 fprintf (stderr,"\n\n wrong interval of indices in function stress (%s, line %d).\n",__FILE__,__LINE__);
01380 abort ();
01381 }
01382
01383 Mt->give_elemnodes (eid,nodes);
01384 k=0;
01385 for (i=fi;i<fi+ncomp;i++){
01386 for (j=0;j<nne;j++){
01387 nodval[j]=Mt->nodes[nodes[j]].stress[lcid*tncomp+i];
01388 }
01389 sig[k]=approx (xi,eta,zeta,nodval);
01390 k++;
01391 }
01392 }
01393
01394
01395
01396 void linhex::stresses (long lcid,long eid,long ri,long ci)
01397 {
01398 long i,naep,ncp,sid;
01399
01400 vector coord,sig;
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416 switch (Mm->stre.tape[eid]){
01417 case nowhere:{
01418 break;
01419 }
01420 case intpts:{
01421
01422 break;
01423 }
01424 case enodes:{
01425 break;
01426 }
01427 case userdefined:{
01428
01429 naep = Mm->stre.give_naep (eid);
01430 ncp = Mm->stre.give_ncomp (eid);
01431 sid = Mm->stre.give_sid (eid);
01432 allocv (ncp,sig);
01433 allocv (3,coord);
01434 for (i=0;i<naep;i++){
01435 Mm->stre.give_aepcoord (sid,i,coord);
01436
01437 if (Mp->stressaver==0)
01438
01439 if (Mp->stressaver==1)
01440 appstress (lcid,eid,coord[0],coord[1],coord[2],0,ncp,sig);
01441
01442 Mm->stre.storevalues(lcid,eid,i,sig);
01443 }
01444 destrv (sig);
01445 destrv (coord);
01446 break;
01447 }
01448 default:{
01449 fprintf (stderr,"\n\n unknown stress point is required in function planeelemlq::stresses (%s, line %d).\n",__FILE__,__LINE__);
01450 }
01451 }
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466 }
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476 void linhex::nod_eqother_ip (long eid,long ri,long ci)
01477 {
01478 long i,j,ipp,ncompo;
01479 ivector ipnum(nne),nod(nne);
01480 vector eqother;
01481
01482
01483
01484 ipp=Mt->elements[eid].ipp[ri][ci];
01485 nodip_linhex (ipp,intordsm[0][0],ipnum);
01486
01487
01488 Mt->give_elemnodes (eid,nod);
01489
01490 for (i=0;i<nne;i++){
01491
01492 ncompo = Mm->givencompeqother (ipnum[i],0);
01493 allocv (ncompo,eqother);
01494
01495 Mm->giveeqother (ipnum[i],0,ncompo,eqother.a);
01496
01497
01498 j=nod[i];
01499 Mt->nodes[j].storeother (0,ncompo,eqother);
01500
01501 destrv (eqother);
01502 }
01503 }
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530 void linhex::gl_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x, vector &y, vector &z)
01531 {
01532 integratedquant iq;
01533 iq=locstress;
01534
01535
01536 compute_nlstress (lcid,eid,ri,ci);
01537
01538
01539 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01540 }
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553 void linhex::gnl_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
01554 {
01555 long i,j,k,ipp;
01556 double xi,eta,zeta,jac;
01557 vector w,gp,x(nne),y(nne),z(nne),sig(tncomp),contr(ndofe),r(ndofe);
01558 matrix gm(tncomp,ndofe);
01559
01560
01561 Mt->give_node_coord3d (x,y,z,eid);
01562
01563 eldispl (lcid,eid,r.a);
01564
01565
01566 fillv (0.0,ifor);
01567
01568
01569 allocv (intordsm[0][0],gp);
01570
01571 allocv (intordsm[0][0],w);
01572
01573
01574 gauss_points (gp.a,w.a,intordsm[0][0]);
01575
01576
01577 ipp=Mt->elements[eid].ipp[ri][ci];
01578
01579
01580 for (i=0;i<intordsm[0][0];i++){
01581 xi=gp[i];
01582 for (j=0;j<intordsm[0][0];j++){
01583 eta=gp[j];
01584 for (k=0;k<intordsm[0][0];k++){
01585 zeta=gp[k];
01586
01587
01588 if (Mp->strcomp==1)
01589 Mm->computenlstresses (ipp);
01590
01591 if (Mp->nodeintfor==1){
01592 Mm->givestress (lcid,ipp,sig);
01593
01594
01595 gngeom_matrix (gm,r,x,y,z,xi,eta,zeta,jac);
01596
01597 mtxv (gm,sig,contr);
01598
01599 cmulv (jac*w[i]*w[j]*w[k],contr);
01600
01601 addv (contr,ifor,ifor);
01602 }
01603
01604 ipp++;
01605 }
01606 }
01607 }
01608 destrv (w); destrv (gp);
01609 }
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624 void linhex::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y,vector &z)
01625 {
01626 integratedquant iq;
01627 iq=nonlocstress;
01628
01629
01630 compute_nonloc_nlstress (lcid,eid,ri,ci);
01631
01632
01633 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01634
01635 }
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650 void linhex::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y,vector &z)
01651 {
01652 integratedquant iq;
01653 iq=stressincr;
01654
01655
01656 compute_nlstressincr (lcid,eid,ri,ci);
01657
01658 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
01659 }
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674 void linhex::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y,vector &z)
01675 {
01676 integratedquant iq;
01677 iq=eigstress;
01678
01679
01680 if (Mp->eigstrcomp)
01681 compute_eigstress (lcid,eid,ri,ci);
01682
01683
01684 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y,z);
01685 }
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698 void linhex::res_internal_forces (long lcid,long eid,vector &ifor)
01699 {
01700 long transf;
01701 ivector nodes (nne);
01702 vector v(ndofe),x(nne),y(nne),z(nne);
01703
01704 Mt->give_node_coord3d (x,y,z,eid);
01705
01706 gl_internal_forces (lcid,eid,0,0,ifor,x,y,z);
01707
01708
01709
01710
01711
01712 Mt->give_elemnodes (eid,nodes);
01713 transf = Mt->locsystems (nodes);
01714 if (transf>0){
01715 matrix tmat (ndofe,ndofe);
01716 transf_matrix (nodes,tmat);
01717
01718 glvectortransf (ifor,v,tmat);
01719 copyv (v,ifor);
01720 }
01721 }
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734 void linhex::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
01735 {
01736 long transf;
01737 ivector nodes (nne);
01738 vector v(ndofe),x(nne),y(nne),z(nne);
01739
01740 Mt->give_node_coord3d (x,y,z,eid);
01741
01742 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y,z);
01743
01744
01745
01746 Mt->give_elemnodes (eid,nodes);
01747 transf = Mt->locsystems (nodes);
01748 if (transf>0){
01749 matrix tmat (ndofe,ndofe);
01750 transf_matrix (nodes,tmat);
01751
01752 glvectortransf (ifor,v,tmat);
01753 copyv (v,ifor);
01754 }
01755 }
01756
01757
01758
01759 void linhex::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01760 {
01761 long transf;
01762 ivector nodes (nne);
01763 vector v(ndofe),x(nne),y(nne),z(nne);
01764
01765 Mt->give_node_coord3d (x,y,z,eid);
01766
01767 incr_internal_forces (lcid,eid,0,0,ifor,x,y,z);
01768
01769
01770
01771 Mt->give_elemnodes (eid,nodes);
01772 transf = Mt->locsystems (nodes);
01773 if (transf>0){
01774 matrix tmat (ndofe,ndofe);
01775 transf_matrix (nodes,tmat);
01776
01777 glvectortransf (ifor,v,tmat);
01778 copyv (v,ifor);
01779 }
01780 }
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792 void linhex::res_eigstrain_forces (long lcid,long eid,vector &nfor)
01793 {
01794 long transf;
01795 ivector nodes (nne);
01796 vector v(ndofe),x(nne),y(nne),z(nne);
01797
01798 Mt->give_node_coord3d (x,y,z,eid);
01799
01800 eigstrain_forces (lcid,eid,0,0,nfor,x,y,z);
01801
01802
01803
01804 Mt->give_elemnodes (eid,nodes);
01805 transf = Mt->locsystems (nodes);
01806 if (transf>0){
01807 matrix tmat (ndofe,ndofe);
01808 transf_matrix (nodes,tmat);
01809
01810 glvectortransf (nfor,v,tmat);
01811 copyv (v,nfor);
01812 }
01813 }
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826 void linhex::compute_nlstress (long lcid,long eid,long ri,long ci)
01827 {
01828 long i,j,k,ipp;
01829
01830 ipp=Mt->elements[eid].ipp[ri][ci];
01831
01832 for (i=0;i<intordsm[0][0];i++){
01833 for (j=0;j<intordsm[0][0];j++){
01834 for (k=0;k<intordsm[0][0];k++){
01835
01836
01837 if (Mp->strcomp==1)
01838 Mm->computenlstresses (ipp);
01839
01840 ipp++;
01841 }
01842 }
01843 }
01844 }
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857 void linhex::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01858 {
01859 long i,j,k,ipp;
01860
01861 ipp=Mt->elements[eid].ipp[ri][ci];
01862
01863 for (i=0;i<intordsm[0][0];i++){
01864 for (j=0;j<intordsm[0][0];j++){
01865 for (k=0;k<intordsm[0][0];k++){
01866
01867
01868 if (Mp->strcomp==1)
01869 Mm->computenlstressesincr (ipp);
01870
01871 ipp++;
01872 }
01873 }
01874 }
01875 }
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889 void linhex::local_values (long lcid,long eid,long ri,long ci)
01890 {
01891 long i,j,k,ipp;
01892
01893 ipp=Mt->elements[eid].ipp[ri][ci];
01894
01895 for (i=0;i<intordsm[0][0];i++){
01896 for (j=0;j<intordsm[0][0];j++){
01897 for (k=0;k<intordsm[0][0];k++){
01898
01899
01900 if (Mp->strcomp==1)
01901 Mm->computenlstresses (ipp);
01902
01903 ipp++;
01904 }
01905 }
01906 }
01907 }
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920 void linhex::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01921 {
01922 long i,j,k,ipp;
01923
01924 ipp=Mt->elements[eid].ipp[ri][ci];
01925 for (i=0;i<intordsm[0][0];i++){
01926 for (j=0;j<intordsm[0][0];j++){
01927 for (k=0;k<intordsm[0][0];k++){
01928
01929 if (Mp->strcomp==1)
01930 Mm->compnonloc_nlstresses (ipp);
01931 ipp++;
01932 }
01933 }
01934 }
01935 }
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948 void linhex::compute_eigstress (long lcid,long eid,long ri,long ci)
01949 {
01950 long i,j,k,ipp;
01951 vector eigstr(tncomp),sig(tncomp);
01952 matrix d(tncomp,tncomp);
01953
01954 ipp=Mt->elements[eid].ipp[ri][ci];
01955
01956 for (i=0;i<intordsm[0][0];i++){
01957 for (j=0;j<intordsm[0][0];j++){
01958 for (k=0;k<intordsm[0][0];k++){
01959 Mm->giveeigstrain (ipp,eigstr);
01960
01961 Mm->matstiff (d,ipp);
01962 mxv (d,eigstr,sig);
01963 Mm->storeeigstress (ipp,sig);
01964 ipp++;
01965 }
01966 }
01967 }
01968 }
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986 void linhex::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x, vector &y, vector &z)
01987 {
01988 long i,j,k,ipp;
01989 double xi,eta,zeta,jac;
01990 vector w,gp,ipv(tncomp),contr(ndofe);
01991 matrix gm(tncomp,ndofe);
01992
01993 fillv (0.0,nv);
01994
01995 allocv (intordsm[0][0],gp);
01996 allocv (intordsm[0][0],w);
01997
01998 gauss_points (gp.a,w.a,intordsm[0][0]);
01999 ipp=Mt->elements[eid].ipp[ri][ci];
02000
02001 for (i=0;i<intordsm[0][0];i++){
02002 xi=gp[i];
02003 for (j=0;j<intordsm[0][0];j++){
02004 eta=gp[j];
02005 for (k=0;k<intordsm[0][0];k++){
02006 zeta=gp[k];
02007
02008 Mm->givequantity (iq,lcid,ipp,cncomp[0],ipv);
02009
02010 geom_matrix (gm,x,y,z,xi,eta,zeta,jac);
02011
02012 mtxv (gm,ipv,contr);
02013 cmulv (jac*w[i]*w[j]*w[k],contr);
02014
02015 addv (contr,nv,nv);
02016 ipp++;
02017 }
02018 }
02019 }
02020 destrv (w); destrv (gp);
02021 }
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040 void linhex::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
02041 {
02042 long i,j,k,ii;
02043 double xi,eta,zeta;
02044 vector x(nne),y(nne),z(nne),w(intordsm[ri][ci]),gp(intordsm[ri][ci]);
02045
02046 gauss_points (gp.a,w.a,intordsm[ri][ci]);
02047 Mt->give_node_coord3d (x,y,z,eid);
02048 ii=Mt->elements[eid].ipp[ri][ci];
02049
02050 for (i=0;i<intordsm[ri][ci];i++){
02051 xi=gp[i];
02052 for (j=0;j<intordsm[ri][ci];j++){
02053 eta=gp[j];
02054 for (k=0;k<intordsm[ri][ci];k++){
02055 zeta=gp[k];
02056 if (ii==ipp){
02057 coord[0]=approx (xi,eta,zeta,x);
02058 coord[1]=approx (xi,eta,zeta,y);
02059 coord[2]=approx (xi,eta,zeta,z);
02060 }
02061 ii++;
02062 }
02063 }
02064 }
02065 }
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091 void linhex::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
02092 {
02093 long i, j, k, l, m, ipp;
02094 long ii, jj, nv = nodval.n;
02095 double xi, eta, zeta, ipval;
02096 vector w, gp, anv(nne);
02097 long nstra, nstre, ncompstr, ncompeqother;
02098 long idstra, idstre, idoth, idic;
02099 inictype ict;
02100
02101 nstra = idstra = nstre = idstre = idoth = idic = 0;
02102
02103 ict = ictn[0];
02104 for (i=0; i<nne; i++)
02105 {
02106 if (ictn[i] != ict)
02107 {
02108 print_err("Incompatible types of initial conditions on element %ld\n"
02109 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
02110 abort();
02111 }
02112 }
02113 for (j = 0; j < nv; j++)
02114 {
02115 for(i = 0; i < nne; i++)
02116 anv[i] = nodval[i][j];
02117 for (ii = 0; ii < nb; ii++)
02118 {
02119 for (jj = 0; jj < nb; jj++)
02120 {
02121 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
02122 if (intordsm[ii][jj] == 0)
02123 continue;
02124 allocv (intordsm[ii][jj],gp);
02125 allocv (intordsm[ii][jj],w);
02126 gauss_points (gp.a,w.a,intordsm[ii][jj]);
02127 for (k = 0; k < intordsm[ii][jj]; k++)
02128 {
02129 xi=gp[k];
02130 for (l = 0; l < intordsm[ii][jj]; l++)
02131 {
02132 eta=gp[l];
02133 for (m = 0; m < intordsm[ii][jj]; m++)
02134 {
02135 zeta=gp[m];
02136
02137 ipval = approx (xi,eta,zeta,anv);
02138 ncompstr = Mm->ip[ipp].ncompstr;
02139 ncompeqother = Mm->ip[ipp].ncompeqother;
02140 if ((ictn[0] & inistrain) && (j < ncompstr))
02141 {
02142 Mm->ip[ipp].strain[idstra] += ipval;
02143 ipp++;
02144 continue;
02145 }
02146 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02147 {
02148 Mm->ip[ipp].stress[idstre] += ipval;
02149 ipp++;
02150 continue;
02151 }
02152 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
02153 {
02154 Mm->ip[ipp].eqother[idoth] += ipval;
02155 ipp++;
02156 continue;
02157 }
02158 if ((ictn[0] & inicond) && (j < nv))
02159 {
02160 if (Mm->ic[ipp] == NULL)
02161 {
02162 Mm->ic[ipp] = new double[nv-j];
02163 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
02164 }
02165 Mm->ic[ipp][idic] += ipval;
02166 ipp++;
02167 continue;
02168 }
02169 ipp++;
02170 }
02171 }
02172 }
02173 destrv (gp); destrv (w);
02174 }
02175 }
02176 ipp=Mt->elements[eid].ipp[ri][ci];
02177 ncompstr = Mm->ip[ipp].ncompstr;
02178 ncompeqother = Mm->ip[ipp].ncompeqother;
02179 if ((ictn[0] & inistrain) && (j < ncompstr))
02180 {
02181 nstra++;
02182 idstra++;
02183 continue;
02184 }
02185 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
02186 {
02187 nstre++;
02188 idstre++;
02189 continue;
02190 }
02191 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
02192 {
02193 idoth++;
02194 continue;
02195 }
02196 if ((ictn[0] & inicond) && (j < nv))
02197 {
02198 idic++;
02199 continue;
02200 }
02201 }
02202 }
02203
02204
02205
02206
02207
02208
02209
02210
02211 void linhex::ipvolume (long eid,long ri,long ci)
02212 {
02213 long i,j,k,ii,jj,ipp;
02214 double xi,eta,zeta,jac;
02215 vector x(nne),y(nne),z(nne),w,gp;
02216
02217 Mt->give_node_coord3d (x,y,z,eid);
02218
02219 for (ii=0;ii<nb;ii++){
02220 for (jj=0;jj<nb;jj++){
02221 if (intordsm[ii][jj]==0) continue;
02222
02223 allocv (intordsm[ii][jj],w);
02224 allocv (intordsm[ii][jj],gp);
02225
02226 gauss_points (gp.a,w.a,intordsm[ii][jj]);
02227
02228 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
02229
02230 for (i=0;i<intordsm[ii][jj];i++){
02231 xi=gp[i];
02232 for (j=0;j<intordsm[ii][jj];j++){
02233 eta=gp[j];
02234 for (k=0;k<intordsm[ii][jj];k++){
02235 zeta=gp[k];
02236
02237 jac_3d (jac,x,y,z,xi,eta,zeta);
02238 jac=fabs(jac);
02239
02240 jac*=w[i]*w[j]*w[k];
02241
02242 Mm->storeipvol (ipp,jac);
02243 ipp++;
02244 }
02245 }
02246 }
02247 destrv (gp); destrv (w);
02248 }
02249 }
02250
02251 }
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263 void linhex::node_forces_surf (long lcid,long eid,long *is,double *nv,vector &nf)
02264 {
02265 long i,j,transf;
02266 double xi,eta,zeta,jac;
02267 double *tnv;
02268 ivector nodes(nne);
02269 vector x(nne),y(nne),z(nne),gp,w,av(ndofe),v(ndofe);
02270 matrix n(napfun,ndofe),an(napfun,ndofe),am(ndofe,ndofe),tmat;
02271
02272 tnv = new double [12];
02273
02274
02275 Mt->give_elemnodes(eid, nodes);
02276
02277 Mt->give_node_coord3d (x,y,z,eid);
02278 allocv (intordb,w);
02279 allocv (intordb,gp);
02280 gauss_points (gp.a,w.a,intordb);
02281
02282
02283
02284 if (is[0]>0 ){
02285 fillm (0.0,am);
02286 xi=1.0;
02287 for (i=0;i<intordb;i++){
02288 eta=gp[i];
02289 for (j=0;j<intordb;j++){
02290 zeta=gp[j];
02291
02292 jac2d_3d (jac,x,y,z,eta,zeta,0);
02293 bf_matrix (n,xi,eta,zeta);
02294 jac = jac*w[i]*w[j];
02295 nnj (am.a,n.a,jac,n.m,n.n);
02296 }
02297 }
02298
02299 if (is[0]==1){
02300 av[0] = nv[0*3+0];
02301 av[1] = nv[0*3+1];
02302 av[2] = nv[0*3+2];
02303
02304 av[9] = nv[1*3+0];
02305 av[10] = nv[1*3+1];
02306 av[11] = nv[1*3+2];
02307
02308 av[21] = nv[2*3+0];
02309 av[22] = nv[2*3+1];
02310 av[23] = nv[2*3+2];
02311
02312 av[12] = nv[3*3+0];
02313 av[13] = nv[3*3+1];
02314 av[14] = nv[3*3+2];
02315 }
02316
02317
02318 if (is[0]==2){
02319
02320 av[0] = nv[0*3+0];
02321 av[1] = nv[0*3+1];
02322 av[2] = nv[0*3+2];
02323
02324 av[3] = nv[1*3+0];
02325 av[4] = nv[1*3+1];
02326 av[5] = nv[1*3+2];
02327
02328 av[6] = nv[2*3+0];
02329 av[7] = nv[2*3+1];
02330 av[8] = nv[2*3+2];
02331
02332 av[9] = nv[3*3+0];
02333 av[10] = nv[3*3+1];
02334 av[11] = nv[3*3+2];
02335
02336 locglob_nodeval (0,av,tnv,x,y,z);
02337 fillv(0.0,av);
02338
02339 av[0] = tnv[0*3+0];
02340 av[1] = tnv[0*3+1];
02341 av[2] = tnv[0*3+2];
02342
02343 av[9] = tnv[1*3+0];
02344 av[10] = tnv[1*3+1];
02345 av[11] = tnv[1*3+2];
02346
02347 av[21] = tnv[2*3+0];
02348 av[22] = tnv[2*3+1];
02349 av[23] = tnv[2*3+2];
02350
02351 av[12] = tnv[3*3+0];
02352 av[13] = tnv[3*3+1];
02353 av[14] = tnv[3*3+2];
02354 }
02355
02356
02357 mxv (am,av,v);
02358 addv (v,nf,nf);
02359 }
02360
02361
02362
02363
02364 if (is[1]>0 ){
02365 fillm (0.0,am);
02366 eta=1.0;
02367 for (i=0;i<intordb;i++){
02368 xi=gp[i];
02369 for (j=0;j<intordb;j++){
02370 zeta=gp[j];
02371
02372 jac2d_3d (jac,x,y,z,xi,zeta,1);
02373 bf_matrix (n,xi,eta,zeta);
02374 jac = jac*w[i]*w[j];
02375 nnj (am.a,n.a,jac,n.m,n.n);
02376 }
02377 }
02378
02379 if (is[1]==1){
02380 av[3] = nv[12+0*3+0];
02381 av[4] = nv[12+0*3+1];
02382 av[5] = nv[12+0*3+2];
02383
02384 av[0] = nv[12+1*3+0];
02385 av[1] = nv[12+1*3+1];
02386 av[2] = nv[12+1*3+2];
02387
02388 av[12] = nv[12+2*3+0];
02389 av[13] = nv[12+2*3+1];
02390 av[14] = nv[12+2*3+2];
02391
02392 av[15] = nv[12+3*3+0];
02393 av[16] = nv[12+3*3+1];
02394 av[17] = nv[12+3*3+2];
02395 }
02396
02397
02398 if (is[1]==2){
02399
02400 av[0] = nv[12+0*3+0];
02401 av[1] = nv[12+0*3+1];
02402 av[2] = nv[12+0*3+2];
02403
02404 av[3] = nv[12+1*3+0];
02405 av[4] = nv[12+1*3+1];
02406 av[5] = nv[12+1*3+2];
02407
02408 av[6] = nv[12+2*3+0];
02409 av[7] = nv[12+2*3+1];
02410 av[8] = nv[12+2*3+2];
02411
02412 av[9] = nv[12+3*3+0];
02413 av[10] = nv[12+3*3+1];
02414 av[11] = nv[12+3*3+2];
02415
02416 locglob_nodeval (1,av,tnv,x,y,z);
02417 fillv(0.0,av);
02418
02419 av[3] = tnv[0*3+0];
02420 av[4] = tnv[0*3+1];
02421 av[5] = tnv[0*3+2];
02422
02423 av[0] = tnv[1*3+0];
02424 av[1] = tnv[1*3+1];
02425 av[2] = tnv[1*3+2];
02426
02427 av[12] = tnv[2*3+0];
02428 av[13] = tnv[2*3+1];
02429 av[14] = tnv[2*3+2];
02430
02431 av[15] = tnv[3*3+0];
02432 av[16] = tnv[3*3+1];
02433 av[17] = tnv[3*3+2];
02434
02435 }
02436
02437 mxv (am,av,v);
02438 addv (v,nf,nf);
02439 }
02440
02441
02442
02443 if (is[2]>0 ){
02444 fillm (0.0,am);
02445 xi=-1.0;
02446 for (i=0;i<intordb;i++){
02447 eta=gp[i];
02448 for (j=0;j<intordb;j++){
02449 zeta=gp[j];
02450
02451 jac2d_3d (jac,x,y,z,eta,zeta,2);
02452 bf_matrix (n,xi,eta,zeta);
02453 jac = jac*w[i]*w[j];
02454 nnj (am.a,n.a,jac,n.m,n.n);
02455 }
02456 }
02457
02458 if (is[2]==1){
02459 av[6] = nv[24+0*3+0];
02460 av[7] = nv[24+0*3+1];
02461 av[8] = nv[24+0*3+2];
02462
02463 av[3] = nv[24+1*3+0];
02464 av[4] = nv[24+1*3+1];
02465 av[5] = nv[24+1*3+2];
02466
02467 av[15] = nv[24+2*3+0];
02468 av[16] = nv[24+2*3+1];
02469 av[17] = nv[24+2*3+2];
02470
02471 av[18] = nv[24+3*3+0];
02472 av[19] = nv[24+3*3+1];
02473 av[20] = nv[24+3*3+2];
02474 }
02475
02476
02477 if (is[2]==2){
02478
02479 av[0] = nv[24+0*3+0];
02480 av[1] = nv[24+0*3+1];
02481 av[2] = nv[24+0*3+2];
02482
02483 av[3] = nv[24+1*3+0];
02484 av[4] = nv[24+1*3+1];
02485 av[5] = nv[24+1*3+2];
02486
02487 av[6] = nv[24+2*3+0];
02488 av[7] = nv[24+2*3+1];
02489 av[8] = nv[24+2*3+2];
02490
02491 av[9] = nv[24+3*3+0];
02492 av[10] = nv[24+3*3+1];
02493 av[11] = nv[24+3*3+2];
02494
02495 locglob_nodeval (2,av,tnv,x,y,z);
02496 fillv(0.0,av);
02497
02498 av[6] = tnv[0*3+0];
02499 av[7] = tnv[0*3+1];
02500 av[8] = tnv[0*3+2];
02501
02502 av[3] = tnv[1*3+0];
02503 av[4] = tnv[1*3+1];
02504 av[5] = tnv[1*3+2];
02505
02506 av[15] = tnv[2*3+0];
02507 av[16] = tnv[2*3+1];
02508 av[17] = tnv[2*3+2];
02509
02510 av[18] = tnv[3*3+0];
02511 av[19] = tnv[3*3+1];
02512 av[20] = tnv[3*3+2];
02513
02514 }
02515
02516 mxv (am,av,v);
02517 addv (v,nf,nf);
02518 }
02519
02520
02521
02522
02523 if (is[3]>0 ){
02524 fillm (0.0,am);
02525 eta=-1.0;
02526 for (i=0;i<intordb;i++){
02527 xi=gp[i];
02528 for (j=0;j<intordb;j++){
02529 zeta=gp[j];
02530
02531 jac2d_3d (jac,x,y,z,xi,zeta,3);
02532 bf_matrix (n,xi,eta,zeta);
02533 jac = jac*w[i]*w[j];
02534 nnj (am.a,n.a,jac,n.m,n.n);
02535 }
02536 }
02537
02538 if (is[3]==1){
02539 av[9] = nv[36+0*3+0];
02540 av[10] = nv[36+0*3+1];
02541 av[11] = nv[36+0*3+2];
02542
02543 av[6] = nv[36+1*3+0];
02544 av[7] = nv[36+1*3+1];
02545 av[8] = nv[36+1*3+2];
02546
02547 av[18] = nv[36+2*3+0];
02548 av[19] = nv[36+2*3+1];
02549 av[20] = nv[36+2*3+2];
02550
02551 av[21] = nv[36+3*3+0];
02552 av[22] = nv[36+3*3+1];
02553 av[23] = nv[36+3*3+2];
02554 }
02555
02556
02557 if (is[3]==2){
02558
02559 av[0] = nv[36+0*3+0];
02560 av[1] = nv[36+0*3+1];
02561 av[2] = nv[36+0*3+2];
02562
02563 av[3] = nv[36+1*3+0];
02564 av[4] = nv[36+1*3+1];
02565 av[5] = nv[36+1*3+2];
02566
02567 av[6] = nv[36+2*3+0];
02568 av[7] = nv[36+2*3+1];
02569 av[8] = nv[36+2*3+2];
02570
02571 av[9] = nv[36+3*3+0];
02572 av[10] = nv[36+3*3+1];
02573 av[11] = nv[36+3*3+2];
02574
02575 locglob_nodeval (3,av,tnv,x,y,z);
02576 fillv(0.0,av);
02577
02578 av[9] = tnv[0*3+0];
02579 av[10]= tnv[0*3+1];
02580 av[11]= tnv[0*3+2];
02581
02582 av[6] = tnv[1*3+0];
02583 av[7] = tnv[1*3+1];
02584 av[8] = tnv[1*3+2];
02585
02586 av[18] = tnv[2*3+0];
02587 av[19] = tnv[2*3+1];
02588 av[20] = tnv[2*3+2];
02589
02590 av[21] = tnv[3*3+0];
02591 av[22] = tnv[3*3+1];
02592 av[23] = tnv[3*3+2];
02593
02594 }
02595
02596 mxv (am,av,v);
02597 addv (v,nf,nf);
02598 }
02599
02600
02601
02602 if (is[4]>0 ){
02603 fillm (0.0,am);
02604 zeta=1.0;
02605 for (i=0;i<intordb;i++){
02606 xi=gp[i];
02607 for (j=0;j<intordb;j++){
02608 eta=gp[j];
02609
02610 jac2d_3d (jac,x,y,z,xi,eta,4);
02611 bf_matrix (n,xi,eta,zeta);
02612 jac = jac*w[i]*w[j];
02613 nnj (am.a,n.a,jac,n.m,n.n);
02614 }
02615 }
02616
02617 if (is[4]==1){
02618 av[0] = nv[48+0*3+0];
02619 av[1] = nv[48+0*3+1];
02620 av[2] = nv[48+0*3+2];
02621
02622 av[3] = nv[48+1*3+0];
02623 av[4] = nv[48+1*3+1];
02624 av[5] = nv[48+1*3+2];
02625
02626 av[6] = nv[48+2*3+0];
02627 av[7] = nv[48+2*3+1];
02628 av[8] = nv[48+2*3+2];
02629
02630 av[9] = nv[48+3*3+0];
02631 av[10] = nv[48+3*3+1];
02632 av[11] = nv[48+3*3+2];
02633 }
02634
02635
02636 if (is[4]==2){
02637
02638 av[0] = nv[48+0*3+0];
02639 av[1] = nv[48+0*3+1];
02640 av[2] = nv[48+0*3+2];
02641
02642 av[3] = nv[48+1*3+0];
02643 av[4] = nv[48+1*3+1];
02644 av[5] = nv[48+1*3+2];
02645
02646 av[6] = nv[48+2*3+0];
02647 av[7] = nv[48+2*3+1];
02648 av[8] = nv[48+2*3+2];
02649
02650 av[9] = nv[48+3*3+0];
02651 av[10] = nv[48+3*3+1];
02652 av[11] = nv[48+3*3+2];
02653
02654 locglob_nodeval (4,av,tnv,x,y,z);
02655 fillv(0.0,av);
02656
02657 av[0] = tnv[0*3+0];
02658 av[1] = tnv[0*3+1];
02659 av[2] = tnv[0*3+2];
02660
02661 av[3] = tnv[1*3+0];
02662 av[4] = tnv[1*3+1];
02663 av[5] = tnv[1*3+2];
02664
02665 av[6] = tnv[2*3+0];
02666 av[7] = tnv[2*3+1];
02667 av[8] = tnv[2*3+2];
02668
02669 av[9] = tnv[3*3+0];
02670 av[10] = tnv[3*3+1];
02671 av[11] = tnv[3*3+2];
02672
02673 }
02674
02675 mxv (am,av,v);
02676 addv (v,nf,nf);
02677 }
02678
02679
02680
02681
02682 if (is[5]>0 ){
02683 fillm (0.0,am);
02684 zeta=-1.0;
02685 for (i=0;i<intordb;i++){
02686 xi=gp[i];
02687 for (j=0;j<intordb;j++){
02688 eta=gp[j];
02689
02690 jac2d_3d (jac,x,y,z,xi,eta,5);
02691 bf_matrix (n,xi,eta,zeta);
02692 jac = jac*w[i]*w[j];
02693 nnj (am.a,n.a,jac,n.m,n.n);
02694 }
02695 }
02696
02697 if (is[5]==1){
02698 av[12] = nv[60+0*3+0];
02699 av[13] = nv[60+0*3+1];
02700 av[14] = nv[60+0*3+2];
02701
02702 av[21] = nv[60+1*3+0];
02703 av[22] = nv[60+1*3+1];
02704 av[23] = nv[60+1*3+2];
02705
02706 av[18] = nv[60+2*3+0];
02707 av[19] = nv[60+2*3+1];
02708 av[20] = nv[60+2*3+2];
02709
02710 av[15] = nv[60+3*3+0];
02711 av[16] = nv[60+3*3+1];
02712 av[17] = nv[60+3*3+2];
02713 }
02714
02715
02716 if (is[5]==2){
02717
02718 av[0] = nv[60+0*3+0];
02719 av[1] = nv[60+0*3+1];
02720 av[2] = nv[60+0*3+2];
02721
02722 av[3] = nv[60+1*3+0];
02723 av[4] = nv[60+1*3+1];
02724 av[5] = nv[60+1*3+2];
02725
02726 av[6] = nv[60+2*3+0];
02727 av[7] = nv[60+2*3+1];
02728 av[8] = nv[60+2*3+2];
02729
02730 av[9] = nv[60+3*3+0];
02731 av[10] = nv[60+3*3+1];
02732 av[11] = nv[60+3*3+2];
02733
02734 locglob_nodeval (5,av,tnv,x,y,z);
02735 fillv(0.0,av);
02736
02737 av[12] = tnv[0*3+0];
02738 av[13] = tnv[0*3+1];
02739 av[14] = tnv[0*3+2];
02740
02741 av[21] = tnv[1*3+0];
02742 av[22] = tnv[1*3+1];
02743 av[23] = tnv[1*3+2];
02744
02745 av[18] = tnv[2*3+0];
02746 av[19] = tnv[2*3+1];
02747 av[20] = tnv[2*3+2];
02748
02749 av[15] = tnv[3*3+0];
02750 av[16] = tnv[3*3+1];
02751 av[17] = tnv[3*3+2];
02752
02753 }
02754
02755 mxv (am,av,v);
02756 addv (v,nf,nf);
02757 }
02758
02759 transf = Mt->locsystems (nodes);
02760 if (transf>0)
02761 {
02762 allocm(ndofe,ndofe,tmat);
02763 transf_matrix (nodes,tmat);
02764 lgvectortransfblock (nf,tmat);
02765 }
02766 destrv (gp); destrv (w);
02767 delete [] tnv;
02768 }
02769
02770 void linhex::locglob_nodeval (long is,vector &nv,double *tnv,vector &x,vector &y,vector &z)
02771 {
02772 double norm;
02773 vector g1(3),g2(3),g3(3),lv(3),gv(3);
02774 matrix t(3,3);
02775
02776 if (is==0){
02777 g1[0]=x[4]-x[7];
02778 g1[1]=y[4]-y[7];
02779 g1[2]=z[4]-z[7];
02780
02781 g2[0]=x[3]-x[7];
02782 g2[1]=y[3]-y[7];
02783 g2[2]=z[3]-z[7];
02784 }
02785 if (is==1){
02786 g1[0]=x[5]-x[4];
02787 g1[1]=y[5]-y[4];
02788 g1[2]=z[5]-z[4];
02789
02790 g2[0]=x[0]-x[4];
02791 g2[1]=y[0]-y[4];
02792 g2[2]=z[0]-z[4];
02793 }
02794 if (is==2){
02795 g1[0]=x[6]-x[5];
02796 g1[1]=y[6]-y[5];
02797 g1[2]=z[6]-z[5];
02798
02799 g2[0]=x[1]-x[5];
02800 g2[1]=y[1]-y[5];
02801 g2[2]=z[1]-z[5];
02802 }
02803 if (is==3){
02804 g1[0]=x[7]-x[6];
02805 g1[1]=y[7]-y[6];
02806 g1[2]=z[7]-z[6];
02807
02808 g2[0]=x[2]-x[6];
02809 g2[1]=y[2]-y[6];
02810 g2[2]=z[2]-z[6];
02811 }
02812 if (is==4){
02813 g1[0]=x[3]-x[2];
02814 g1[1]=y[3]-y[2];
02815 g1[2]=z[3]-z[2];
02816
02817 g2[0]=x[1]-x[2];
02818 g2[1]=y[1]-y[2];
02819 g2[2]=z[1]-z[2];
02820 }
02821 if (is==5){
02822 g1[0]=x[5]-x[6];
02823 g1[1]=y[5]-y[6];
02824 g1[2]=z[5]-z[6];
02825
02826 g2[0]=x[7]-x[6];
02827 g2[1]=y[7]-y[6];
02828 g2[2]=z[7]-z[6];
02829 }
02830
02831
02832 scprd (g1,g1,norm);
02833 norm=1.0/sqrt(norm);
02834 cmulv (norm,g1,g1);
02835
02836 scprd (g1,g2,norm);
02837 g2[0]=g2[0]-norm*g1[0];
02838 g2[1]=g2[1]-norm*g1[1];
02839 g2[2]=g2[2]-norm*g1[2];
02840
02841 scprd (g2,g2,norm);
02842 norm=1.0/sqrt(norm);
02843 cmulv (norm,g2,g2);
02844
02845 g3[0]=g1[1]*g2[2]-g1[2]*g2[1];
02846 g3[1]=g1[2]*g2[0]-g1[0]*g2[2];
02847 g3[2]=g1[0]*g2[1]-g1[1]*g2[0];
02848
02849 t[0][0]=g1[0];
02850 t[1][0]=g1[1];
02851 t[2][0]=g1[2];
02852
02853 t[0][1]=g2[0];
02854 t[1][1]=g2[1];
02855 t[2][1]=g2[2];
02856
02857 t[0][2]=g3[0];
02858 t[1][2]=g3[1];
02859 t[2][2]=g3[2];
02860
02861 mxv (t.a,nv.a,tnv,3,3);
02862 mxv (t.a,nv.a+3,tnv+3,3,3);
02863 mxv (t.a,nv.a+6,tnv+6,3,3);
02864 mxv (t.a,nv.a+9,tnv+9,3,3);
02865
02866 }
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878 void linhex::node_forces_surf_old (long lcid,long eid,long *is,double *nv,vector &nf)
02879 {
02880 long i,j,i1,ii;
02881 double xi=0.0,eta=0.0,zeta=0.0,jac, w1,w2;
02882 ivector nodes(nne);
02883 vector gx(nne),gy(nne),gz(nne),x(nnsurf),y(nnsurf),z(nnsurf),gp(intordb),w(intordb),av(ndofe),v(ndofe),pom(3);
02884 matrix n(napfun,ndofe),an(napfun,ndofe),am(ndofe,ndofe), tran(3,3);
02885
02886 Mt->give_elemnodes (eid,nodes);
02887 Mt->give_node_coord3d (gx,gy,gz,eid);
02888 gauss_points (gp.a,w.a,intordb);
02889 fillv (0.0,nf);
02890 fillm (0.0,an);
02891 for (ii=0;ii<6;ii++){
02892 fillv (0.0,av);
02893
02894 if (is[ii] ==0){}
02895 else {
02896 tran_mat( tran, gx, gy, gz, ii);
02897
02898 if (ii ==0){
02899 xi=1.0;
02900 x[0]=gx[0];x[1]=gx[3];x[2]=gx[7];x[3]=gx[4];
02901 y[0]=gy[0];y[1]=gy[3];y[2]=gy[7];y[3]=gy[4];
02902 z[0]=gz[0];z[1]=gz[3];z[2]=gz[7];z[3]=gz[4];
02903 for (i=0;i<intordb;i++){
02904 eta=gp[i]; w1=w[i];
02905 for (j=0;j<intordb;j++){
02906 zeta=gp[j]; w2=w[j];
02907 bf_matrix (n,xi,eta,zeta);
02908 jac2d3d (jac,x,y,z,eta,zeta);
02909 jac = jac*w1*w2;
02910 nnj (am.a,n.a,jac,n.m,n.n);
02911 }
02912 }
02913 for (i=0;i<3;i++){
02914 av[i]=nv[i]; av[9+i]=nv[3+i]; av[21+i]=nv[6+i]; av[12+i]=nv[9+i];
02915 }
02916
02917 if (is[ii] ==2){
02918 lgvectortransfblock (av, tran);
02919 }
02920 mxv (am,av,v);
02921 }
02922
02923 else if (ii ==1){
02924 eta=1.0;
02925 x[0]=gx[1];x[1]=gx[0];x[2]=gx[4];x[3]=gx[5];
02926 y[0]=gy[1];y[1]=gy[0];y[2]=gy[4];y[3]=gy[5];
02927 z[0]=gz[1];z[1]=gz[0];z[2]=gz[4];z[3]=gz[5];
02928 for (i=0;i<intordb;i++){
02929 xi=gp[i]; w1=w[i];
02930 for (j=0;j<intordb;j++){
02931 zeta=gp[j]; w2=w[j];
02932 bf_matrix (n,xi,eta,zeta);
02933 jac2d3d (jac,x,y,z,xi,zeta);
02934 jac = jac*w1*w2;
02935 nnj (am.a,n.a,jac,n.m,n.n);
02936 }
02937 }
02938 for (i=0;i<3;i++){
02939 av[3+i]=nv[12+i]; av[i]=nv[15+i]; av[12+i]=nv[18+i]; av[15+i]=nv[21+i];
02940 }
02941
02942 if (is[ii] ==2){
02943 lgvectortransfblock (av, tran);
02944 }
02945 mxv (am,av,v);
02946 }
02947
02948 else if (ii ==2){
02949 xi=-1.0;
02950 x[0]=gx[2];x[1]=gx[1];x[2]=gx[5];x[3]=gx[6];
02951 y[0]=gy[2];y[1]=gy[1];y[2]=gy[5];y[3]=gy[6];
02952 z[0]=gz[2];z[1]=gz[1];z[2]=gz[5];z[3]=gz[6];
02953 for (i=0;i<intordb;i++){
02954 eta=gp[i]; w1=w[i];
02955 for (j=0;j<intordb;j++){
02956 zeta=gp[j]; w2=w[j];
02957 bf_matrix (n,xi,eta,zeta);
02958 jac2d3d (jac,x,y,z,eta,zeta);
02959 jac = jac*w1*w2;
02960 nnj (am.a,n.a,jac,n.m,n.n);
02961 }
02962 }
02963 for (i=0;i<3;i++){
02964 av[6+i]=nv[24+i]; av[3+i]=nv[27+i]; av[15+i]=nv[30+i]; av[18+i]=nv[33+i];
02965 }
02966
02967 if (is[ii] ==2){
02968 lgvectortransfblock (av, tran);
02969 }
02970 mxv (am,av,v);
02971 }
02972
02973 else if (ii ==3){
02974 eta=-1.0;
02975 x[0]=gx[3];x[1]=gx[2];x[2]=gx[6];x[3]=gx[7];
02976 y[0]=gy[3];y[1]=gy[2];y[2]=gy[6];y[3]=gy[7];
02977 z[0]=gz[3];z[1]=gz[2];z[2]=gz[6];z[3]=gz[7];
02978 for (i=0;i<intordb;i++){
02979 xi=gp[i]; w1=w[i];
02980 for (j=0;j<intordb;j++){
02981 zeta=gp[j]; w2=w[j];
02982 bf_matrix (n,xi,eta,zeta);
02983 jac2d3d (jac,x,y,z,xi,zeta);
02984 jac = jac*w1*w2;
02985
02986 nnj (am.a,n.a,jac,n.m,n.n);
02987 }
02988 }
02989 for (i=0;i<3;i++){
02990 av[9+i]=nv[36+i]; av[6+i]=nv[39+i]; av[18+i]=nv[42+i]; av[21+i]=nv[45+i];
02991 }
02992
02993 if (is[ii] ==2){
02994 lgvectortransfblock (av, tran);
02995 }
02996 mxv (am,av,v);
02997 }
02998
02999 else if (ii ==4) {
03000 zeta=1.0;
03001 x[0]=gx[0];x[1]=gx[1];x[2]=gx[2];x[3]=gx[3];
03002 y[0]=gy[0];y[1]=gy[1];y[2]=gy[2];y[3]=gy[3];
03003 z[0]=gz[0];z[1]=gz[1];z[2]=gz[2];z[3]=gz[3];
03004 for (i=0;i<intordb;i++){
03005 xi=gp[i]; w1=w[i];
03006 for (j=0;j<intordb;j++){
03007 eta=gp[j]; w2=w[j];
03008 bf_matrix (n,xi,eta,zeta);
03009 jac2d3d (jac,x,y,z,xi,eta);
03010 jac = jac*w1*w2;
03011
03012
03013
03014 nnj (am.a,n.a,jac,n.m,n.n);
03015 }
03016 }
03017 for (i=0;i<3;i++){
03018 av[i]=nv[48+i]; av[3+i]=nv[51+i]; av[6+i]=nv[54+i]; av[9+i]=nv[57+i];
03019 }
03020
03021
03022
03023
03024
03025
03026
03027
03028 if (is[ii] ==2){
03029 lgvectortransfblock (av, tran);
03030 }
03031 mxv (am,av,v);
03032 }
03033
03034 else if (ii ==5) {
03035 zeta=-1.0;
03036 x[0]=gx[5];x[1]=gx[4];x[2]=gx[7];x[3]=gx[6];
03037 y[0]=gy[5];y[1]=gy[4];y[2]=gy[7];y[3]=gy[6];
03038 z[0]=gz[5];z[1]=gz[4];z[2]=gz[7];z[3]=gz[6];
03039 for (i=0;i<intordb;i++){
03040 xi=gp[i]; w1=w[i];
03041 for (j=0;j<intordb;j++){
03042 eta=gp[j]; w2=w[j];
03043 bf_matrix (n,xi,eta,zeta);
03044 jac2d3d (jac,x,y,z,xi,eta);
03045 jac = jac*w1*w2;
03046 nnj (am.a,n.a,jac,n.m,n.n);
03047 }
03048 }
03049 for (i=0;i<3;i++){
03050 av[12+i]=nv[60+i]; av[15+i]=nv[63+i]; av[18+i]=nv[66+i]; av[21+i]=nv[69+i];
03051 }
03052
03053 if (is[ii] ==2){
03054 lgvectortransfblock (av, tran);
03055 }
03056 mxv (am,av,v);
03057 }
03058
03059 fprintf (Out,"\n\n zatiz");
03060 for (i=0;i<nne;i++){
03061 i1=3*i;
03062 fprintf (Out,"\n %4ld %20.10e %20.10e %20.10e",i,av[i1],av[i1+1],av[i1+2]);}
03063 fprintf (Out,"\n\n vloc");
03064 for (i=0;i<nne;i++){
03065 i1=3*i;
03066 fprintf (Out,"\n %4ld %20.10e %20.10e %20.10e",i,v[i1],v[i1+1],v[i1+2]);}
03067
03068 fprintf (Out,"\n\n vglob");
03069 for (i=0;i<nne;i++){
03070 i1=3*i;
03071 fprintf (Out,"\n %4ld %20.10e %20.10e %20.10e",i,v[i1],v[i1+1],v[i1+2]);}
03072 addv (nf,v,nf);
03073 }
03074 }
03075
03076
03077
03078
03079
03080
03081
03082
03083
03084
03085
03086
03087 }
03088
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099 void linhex::tran_mat(matrix &tran, vector &gx, vector &gy, vector &gz, long is)
03100 {
03101 long i,i1;
03102 double dl;
03103 matrix a(3,3);
03104
03105 if (is ==4) {
03106 for (i=0; i<3; i++) {
03107 i1=i+1; if(i1>2)i1=i1-3;
03108 a[0][i]=gx[i1]-gx[i];
03109 a[1][i]=gy[i1]-gy[i];
03110 a[2][i]=gz[i1]-gz[i];
03111 }
03112 }
03113
03114 else if (is ==5) {
03115 for (i=4; i<7; i++) {
03116 i1=i+1; if(i1>6)i1=i1-3;
03117 a[0][i-4]=gx[i1]-gx[i];
03118 a[1][i-4]=gy[i1]-gy[i];
03119 a[2][i-4]=gz[i1]-gz[i];
03120 }
03121 }
03122
03123
03124
03125
03126 else {
03127 if(is>0)i1=is-4;else i1=is;
03128 a[0][0]=gx[i1+3]-gx[is];
03129 a[1][0]=gy[i1+3]-gy[is];
03130 a[2][0]=gz[i1+3]-gz[is];
03131 a[0][1]=gx[i1+7]-gx[i1+3];
03132 a[1][1]=gy[i1+7]-gy[i1+3];
03133 a[2][1]=gz[i1+7]-gz[i1+3];
03134 a[0][2]=gx[is]-gx[i1+7];
03135 a[1][2]=gy[is]-gy[i1+7];
03136 a[2][2]=gz[is]-gz[i1+7];
03137 }
03138
03139 dl=sqrt(a[0][0]*a[0][0]+a[1][0]*a[1][0]+a[2][0]*a[2][0]);
03140
03141 tran[0][0]=a[0][0]/dl;
03142 tran[1][0]=a[1][0]/dl;
03143 tran[2][0]=a[2][0]/dl;
03144
03145 tran[0][2]=a[1][0]*a[2][1]-a[2][0]*a[1][1];
03146 tran[1][2]=a[2][0]*a[0][1]-a[0][0]*a[2][1];
03147 tran[2][2]=a[0][0]*a[1][1]-a[1][0]*a[0][1];
03148
03149 dl=sqrt(tran[0][2]*tran[0][2]+tran[1][2]*tran[1][2]+tran[2][2]*tran[2][2]);
03150
03151 tran[0][2]=tran[0][2]/dl;
03152 tran[1][2]=tran[1][2]/dl;
03153 tran[2][2]=tran[2][2]/dl;
03154
03155 tran[0][1]=tran[1][2]*tran[2][0]-tran[2][2]*tran[1][0];
03156 tran[1][1]=tran[2][2]*tran[0][0]-tran[0][2]*tran[2][0];
03157 tran[2][1]=tran[0][2]*tran[1][0]-tran[1][2]*tran[0][0];
03158
03159 dl=sqrt(tran[0][1]*tran[0][1]+tran[1][1]*tran[1][1]+tran[2][1]*tran[2][1]);
03160
03161 tran[0][1]=tran[0][1]/dl;
03162 tran[1][1]=tran[1][1]/dl;
03163 tran[2][1]=tran[2][1]/dl;
03164
03165
03166
03167
03168
03169
03170
03171
03172
03173
03174
03175
03176
03177
03178
03179
03180
03181
03182 }
03183
03184
03185
03186
03187
03188
03189
03190
03191
03192 void linhex::intpointval (long eid,vector &nodval,vector &ipval)
03193 {
03194 long i,j,ii,jj,k,l;
03195 double xi,eta,zeta;
03196 vector w,gp;
03197
03198 l=0;
03199 for (ii=0;ii<nb;ii++){
03200 for (jj=0;jj<nb;jj++){
03201 if (intordsm[ii][jj]==0) continue;
03202
03203 allocv (intordsm[ii][jj],w);
03204 allocv (intordsm[ii][jj],gp);
03205
03206 gauss_points (gp.a,w.a,intordsm[ii][jj]);
03207
03208 for (i=0;i<intordsm[ii][jj];i++){
03209 xi=gp[i];
03210 for (j=0;j<intordsm[ii][jj];j++){
03211 eta=gp[j];
03212 for (k=0;k<intordsm[ii][jj];k++){
03213 zeta=gp[k];
03214
03215 ipval[l]=approx (xi,eta,zeta,nodval);
03216 l++;
03217 }
03218 }
03219 }
03220
03221 destrv (w);
03222 destrv (gp);
03223 }
03224 }
03225 }
03226
03227
03228
03229 void linhex::aver_strains (long lcid,long eid,long ri,long ci,vector &averstra,double &volume)
03230 {
03231 long i,j,k,l,ii,ipp;
03232 double xi,eta,zeta,jac;
03233 vector x(nne),y(nne),z(nne),gp,w,eps;
03234
03235 Mt->give_node_coord3d (x,y,z,eid);
03236
03237 for (ii=0;ii<nb;ii++){
03238 if (intordsm[ii][ii]==0) continue;
03239
03240 allocv (intordsm[ii][ii],gp);
03241 allocv (intordsm[ii][ii],w);
03242 allocv (ncomp[ii],eps);
03243
03244 gauss_points (gp.a,w.a,intordsm[ii][ii]);
03245
03246 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
03247 for (i=0;i<intordsm[ii][ii];i++){
03248 xi=gp[i];
03249 for (j=0;j<intordsm[ii][ii];j++){
03250 eta=gp[j];
03251 for (k=0;k<intordsm[ii][ii];k++){
03252 zeta=gp[k];
03253
03254
03255 Mm->givestress (lcid,ipp,eps);
03256
03257 jac_3d (jac,x,y,z,xi,eta,zeta);
03258 jac=fabs(jac);
03259
03260 volume+=w[i]*w[j]*w[k]*jac;
03261
03262 for (l=0;l<averstra.n;l++){
03263 averstra[l]+=eps[l]*w[i]*w[j]*w[k]*jac;
03264 }
03265
03266 ipp++;
03267 }
03268 }
03269 }
03270
03271
03272 destrv (w); destrv (gp); destrv (eps);
03273 }
03274
03275 }
03276
03277
03278
03279
03280
03281
03282
03283
03284
03285
03286
03287
03288
03289
03290 void linhex::surfnodeval (long surf,vector &nodval,double *list)
03291 {
03292 long i,j,k;
03293 ivector surfnod(nnsurf);
03294
03295 fillv (0.0,nodval);
03296 quadhexahedral_surfnod (surfnod.a,surf);
03297
03298 k=0;
03299 for (i=0;i<nnsurf;i++){
03300 for (j=0;j<napfun;j++){
03301 nodval[k]=list[surf*nnsurf*napfun+k];
03302 k++;
03303 }
03304 }
03305 }
03306
03307
03308
03309
03310
03311
03312
03313
03314
03315
03316
03317 void linhex::find_extreme_strains (vector &min,vector &max,long lcid,long eid)
03318 {
03319 long i,j,ipp;
03320 vector sig(tncomp);
03321
03322
03323 ipp=Mt->elements[eid].ipp[0][0];
03324
03325
03326 for (j=0;j<nip[0][0];j++){
03327
03328
03329 Mm->givestress (lcid,ipp,sig);
03330
03331
03332 for (i=0;i<tncomp;i++){
03333 if (min[i]>sig[i])
03334 min[i]=sig[i];
03335 if (max[i]<sig[i])
03336 max[i]=sig[i];
03337 }
03338
03339 ipp++;
03340 }
03341 }
03342
03343
03344
03345
03346
03347
03348
03349
03350
03351
03352
03353 void linhex::find_extreme_stresses (vector &min,vector &max,long lcid,long eid)
03354 {
03355 long i,j,ipp;
03356 vector sig(tncomp);
03357
03358
03359 ipp=Mt->elements[eid].ipp[0][0];
03360
03361
03362 for (j=0;j<nip[0][0];j++){
03363
03364
03365 Mm->givestress (lcid,ipp,sig);
03366
03367
03368 for (i=0;i<tncomp;i++){
03369 if (min[i]>sig[i])
03370 min[i]=sig[i];
03371 if (max[i]<sig[i])
03372 max[i]=sig[i];
03373 }
03374
03375 ipp++;
03376 }
03377 }