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