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