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