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