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