00001 #include "global.h"
00002 #include "cct.h"
00003 #include "gadaptivity.h"
00004 #include "globmat.h"
00005 #include "genfile.h"
00006 #include "element.h"
00007 #include "node.h"
00008 #include "intpoints.h"
00009 #include <stdlib.h>
00010 #include <math.h>
00011
00012 cctelem::cctelem (void)
00013 {
00014 long i,j;
00015
00016
00017 nne=3;
00018
00019 ndofe=9;
00020
00021 tncomp=5;
00022
00023 napfun=3;
00024
00025 ned=3;
00026
00027 nned=2;
00028
00029 intordmm=3;
00030
00031 ssst=plates;
00032
00033
00034 nb=2;
00035
00036
00037 ncomp = new long [nb];
00038 ncomp[0]=3;
00039 ncomp[1]=2;
00040
00041
00042 cncomp = new long [nb];
00043 cncomp[0]=0;
00044 cncomp[1]=3;
00045
00046
00047
00048 nip = new long* [nb];
00049 intordsm = new long* [nb];
00050 for (i=0;i<nb;i++){
00051 nip[i] = new long [nb];
00052 intordsm[i] = new long [nb];
00053 }
00054
00055 nip[0][0]=1; nip[0][1]=0;
00056 nip[1][0]=0; nip[1][1]=1;
00057
00058
00059 tnip=0;
00060 for (i=0;i<nb;i++){
00061 for (j=0;j<nb;j++){
00062 tnip+=nip[i][j];
00063 }
00064 }
00065
00066 intordsm[0][0]=1; intordsm[0][1]=0;
00067 intordsm[1][0]=0; intordsm[1][1]=1;
00068 }
00069
00070 cctelem::~cctelem (void)
00071 {
00072 long i;
00073
00074 for (i=0;i<nb;i++){
00075 delete [] intordsm[i];
00076 delete [] nip[i];
00077 }
00078 delete [] intordsm;
00079 delete [] nip;
00080
00081 delete [] cncomp;
00082 delete [] ncomp;
00083 }
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 double cctelem::approx (vector &areacoord,vector &nodval)
00095 {
00096 double f;
00097 scprd (areacoord,nodval,f);
00098 return f;
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108 double cctelem::approx_nat (double xi,double eta,vector &nodval)
00109 {
00110 double f;
00111 vector areacoord(3);
00112 areacoord[0]=xi;
00113 areacoord[1]=eta;
00114 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00115 scprd (areacoord,nodval,f);
00116 return f;
00117 }
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129 void cctelem::bf_matrix (matrix &n,vector &x,vector &y,vector &areacoord)
00130 {
00131 vector bf(9),sx(3),sy(3);
00132
00133 sx[0]=x[1]-x[0]; sx[1]=x[2]-x[1]; sx[2]=x[0]-x[2];
00134 sy[0]=y[1]-y[0]; sy[1]=y[2]-y[1]; sy[2]=y[0]-y[2];
00135
00136 bf_cct (bf.a,areacoord.a,sx.a,sy.a);
00137
00138 fillm (0.0,n);
00139
00140 n[0][0]=bf[0];
00141 n[0][1]=bf[1];
00142 n[0][2]=bf[2];
00143 n[0][3]=bf[3];
00144 n[0][4]=bf[4];
00145 n[0][5]=bf[5];
00146 n[0][6]=bf[6];
00147 n[0][7]=bf[7];
00148 n[0][8]=bf[8];
00149
00150 n[1][1]=bf[0];
00151 n[1][4]=bf[3];
00152 n[1][7]=bf[6];
00153
00154 n[2][2]=bf[0];
00155 n[2][5]=bf[3];
00156 n[2][8]=bf[6];
00157
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 void cctelem::geom_matrix_block (matrix &gm,long bi,vector &x,vector &y,vector &areacoord)
00170 {
00171 double det;
00172 vector sx(3),sy(3),b(3),c(3),bf(9),dx(9),dy(9);
00173
00174 sx[0]=x[1]-x[0]; sx[1]=x[2]-x[1]; sx[2]=x[0]-x[2];
00175 sy[0]=y[1]-y[0]; sy[1]=y[2]-y[1]; sy[2]=y[0]-y[2];
00176
00177
00178 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00179
00180 plsb (b.a,y.a,det);
00181 plsc (c.a,x.a,det);
00182
00183 bf_cct (bf.a,areacoord.a,sx.a,sy.a);
00184 dx_cct (dx.a,areacoord.a,b.a,sx.a,sy.a);
00185 dy_cct (dy.a,areacoord.a,c.a,sx.a,sy.a);
00186
00187 fillm (0.0,gm);
00188
00189 if (bi==0){
00190 gm[0][2]=dx[0];
00191 gm[0][5]=dx[3];
00192 gm[0][8]=dx[6];
00193
00194 gm[1][1]=0.0-dy[0];
00195 gm[1][4]=0.0-dy[3];
00196 gm[1][7]=0.0-dy[6];
00197
00198 gm[2][1]=0.0-dx[0];
00199 gm[2][2]=dy[0];
00200 gm[2][4]=0.0-dx[3];
00201 gm[2][5]=dy[3];
00202 gm[2][7]=0.0-dx[6];
00203 gm[2][8]=dy[6];
00204 }
00205 if (bi==1){
00206 gm[0][0]=dy[0];
00207 gm[0][1]=dy[1]-bf[0];
00208 gm[0][2]=dy[2];
00209 gm[0][3]=dy[3];
00210 gm[0][4]=dy[4]-bf[3];
00211 gm[0][5]=dy[5];
00212 gm[0][6]=dy[6];
00213 gm[0][7]=dy[7]-bf[6];
00214 gm[0][8]=dy[8];
00215
00216 gm[1][0]=dx[0];
00217 gm[1][1]=dx[1];
00218 gm[1][2]=dx[2]+bf[0];
00219 gm[1][3]=dx[3];
00220 gm[1][4]=dx[4];
00221 gm[1][5]=dx[5]+bf[3];
00222 gm[1][6]=dx[6];
00223 gm[1][7]=dx[7];
00224 gm[1][8]=dx[8]+bf[6];
00225 }
00226
00227 }
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 void cctelem::geom_matrix (matrix &gm,vector &x,vector &y,vector &areacoord)
00239 {
00240 double det;
00241 vector sx(3),sy(3),b(3),c(3),bf(9),dx(9),dy(9);
00242
00243 sx[0]=x[1]-x[0]; sx[1]=x[2]-x[1]; sx[2]=x[0]-x[2];
00244 sy[0]=y[1]-y[0]; sy[1]=y[2]-y[1]; sy[2]=y[0]-y[2];
00245
00246
00247 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00248
00249 plsb (b.a,y.a,det);
00250 plsc (c.a,x.a,det);
00251
00252 bf_cct (bf.a,areacoord.a,sx.a,sy.a);
00253 dx_cct (dx.a,areacoord.a,b.a,sx.a,sy.a);
00254 dy_cct (dy.a,areacoord.a,c.a,sx.a,sy.a);
00255
00256 fillm (0.0,gm);
00257
00258 gm[0][2]=dx[0];
00259 gm[0][5]=dx[3];
00260 gm[0][8]=dx[6];
00261
00262 gm[1][1]=0.0-dy[0];
00263 gm[1][4]=0.0-dy[3];
00264 gm[1][7]=0.0-dy[6];
00265
00266 gm[2][1]=0.0-dx[0];
00267 gm[2][2]=dy[0];
00268 gm[2][4]=0.0-dx[3];
00269 gm[2][5]=dy[3];
00270 gm[2][7]=0.0-dx[6];
00271 gm[2][8]=dy[6];
00272
00273 gm[3][0]=dy[0];
00274 gm[3][1]=dy[1]-bf[0];
00275 gm[3][2]=dy[2];
00276 gm[3][3]=dy[3];
00277 gm[3][4]=dy[4]-bf[3];
00278 gm[3][5]=dy[5];
00279 gm[3][6]=dy[6];
00280 gm[3][7]=dy[7]-bf[6];
00281 gm[3][8]=dy[8];
00282
00283 gm[4][0]=dx[0];
00284 gm[4][1]=dx[1];
00285 gm[4][2]=dx[2]+bf[0];
00286 gm[4][3]=dx[3];
00287 gm[4][4]=dx[4];
00288 gm[4][5]=dx[5]+bf[3];
00289 gm[4][6]=dx[6];
00290 gm[4][7]=dx[7];
00291 gm[4][8]=dx[8]+bf[6];
00292 }
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 void cctelem::dmatblock (long ri,long ci,matrix &d,matrix &dd,double t)
00305 {
00306 double c;
00307
00308 c=t*t*t;
00309 fillm (0.0,dd);
00310
00311 if (ri==0 && ci==0){
00312 dd[0][0] = c*d[0][0]; dd[0][1] = c*d[0][1]; dd[0][2] = c*d[0][2];
00313 dd[1][0] = c*d[1][0]; dd[1][1] = c*d[1][1]; dd[1][2] = c*d[1][2];
00314 dd[2][0] = c*d[2][0]; dd[2][1] = c*d[2][1]; dd[2][2] = c*d[2][2];
00315 }
00316 if (ri==1 && ci==1){
00317 dd[0][0] = t*d[3][3]*5.0/6.0; dd[0][1] = 0.0;
00318 dd[1][0] = 0.0; dd[1][1] = t*d[4][4]*5.0/6.0;
00319 }
00320
00321 }
00322
00323 void cctelem::dmat (matrix &d,double t)
00324 {
00325 double c;
00326
00327 c=t*t*t;
00328
00329 d[0][0] = c*d[0][0]; d[0][1] = c*d[0][1]; d[0][2] = c*d[0][2];
00330 d[1][0] = c*d[1][0]; d[1][1] = c*d[1][1]; d[1][2] = c*d[1][2];
00331 d[2][0] = c*d[2][0]; d[2][1] = c*d[2][1]; d[2][2] = c*d[2][2];
00332
00333 d[3][3] = t*d[3][3]*5.0/6.0;
00334 d[4][4] = t*d[4][4]*5.0/6.0;
00335 }
00336
00337
00338
00339
00340
00341
00342 void cctelem::transf_matrix (ivector &nodes,matrix &tmat)
00343 {
00344 long i,n,m;
00345 n=nodes.n;
00346 m=tmat.m;
00347 for (i=0;i<m;i++){
00348 tmat[i][i]=1.0;
00349 }
00350
00351 for (i=0;i<n;i++){
00352 if (Mt->nodes[nodes[i]].transf>0){
00353 tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e1[0]; tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e2[0];
00354 tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e1[1]; tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e2[1];
00355 }
00356 }
00357 }
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369 void cctelem::stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
00370 {
00371 long i,ii,jj,ipp;
00372 double jac,det,thick;
00373 ivector nodes(nne);
00374 vector l(3),t(nne),gp1,gp2,w,areacoord(3);
00375 matrix gmr,gmc,dd,d(tncomp,tncomp);
00376
00377 Mt->give_elemnodes (eid,nodes);
00378 Mc->give_thickness (eid,nodes,t);
00379
00380
00381 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00382
00383 fillm (0.0,sm);
00384
00385
00386 for (ii=0;ii<nb;ii++){
00387 allocm (ncomp[ii],ndofe,gmr);
00388 for (jj=0;jj<nb;jj++){
00389 if (intordsm[ii][jj]==0) continue;
00390
00391 allocv (intordsm[ii][jj],w);
00392 allocv (intordsm[ii][jj],gp1);
00393 allocv (intordsm[ii][jj],gp2);
00394
00395 allocm (ncomp[jj],ndofe,gmc);
00396 allocm (ncomp[ii],ncomp[jj],dd);
00397
00398 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
00399
00400 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00401
00402 for (i=0;i<intordsm[ii][jj];i++){
00403 areacoord[0]=gp1[i];
00404 areacoord[1]=gp2[i];
00405 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00406
00407 thick = approx (areacoord,t);
00408
00409
00410 geom_matrix_block (gmr,ii,x,y,areacoord);
00411 geom_matrix_block (gmc,jj,x,y,areacoord);
00412
00413
00414 Mm->matstiff (d,ipp);
00415 dmatblock (ii,jj,d,dd,thick);
00416
00417 jac=w[i]*det;
00418
00419
00420 bdbjac (sm,gmr,dd,gmc,jac);
00421
00422 ipp++;
00423 }
00424
00425 destrm (dd); destrm (gmc); destrv (gp2); destrv (gp1); destrv (w);
00426 }
00427 destrm (gmr);
00428 }
00429
00430 }
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440 void cctelem::res_stiffness_matrix (long eid,matrix &sm)
00441 {
00442 long transf;
00443 ivector nodes(nne);
00444 vector x(nne),y(nne);
00445 Mt->give_node_coord2d (x,y,eid);
00446 Mt->give_elemnodes (eid,nodes);
00447
00448 stiffness_matrix (eid,0,0,sm,x,y);
00449
00450
00451 transf = Mt->locsystems (nodes);
00452 if (transf>0){
00453 matrix tmat (ndofe,ndofe);
00454 transf_matrix (nodes,tmat);
00455 glmatrixtransf (sm,tmat);
00456 }
00457 }
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468 void cctelem::mass_matrix (long eid,matrix &mm,vector &x,vector &y)
00469 {
00470 long i;
00471 double det,ww,thick,rho;
00472 ivector nodes(nne);
00473 vector l(3),gp1(intordmm),gp2(intordmm),w(intordmm),t(nne),dens(nne);
00474 matrix n(3,ndofe);
00475
00476 Mt->give_elemnodes (eid,nodes);
00477 Mc->give_thickness (eid,nodes,t);
00478 Mc->give_density (eid,nodes,dens);
00479
00480 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00481
00482 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00483
00484 fillm (0.0,mm);
00485
00486 for (i=0;i<intordmm;i++){
00487 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00488 ww=w[i];
00489
00490 bf_matrix (n,x,y,l);
00491
00492 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
00493 rho = dens[0]*l[0]+dens[1]*l[1]+dens[2]*l[2];
00494
00495 ww*=det*thick*rho;
00496
00497 nnj (mm.a,n.a,ww,3,ndofe);
00498 }
00499 }
00500
00501 void cctelem::res_mass_matrix (long eid,matrix &mm)
00502 {
00503 vector x(nne),y(nne);
00504 Mt->give_node_coord2d (x,y,eid);
00505 mass_matrix (eid,mm,x,y);
00506 }
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517 void cctelem::load_matrix (long eid,matrix &lm)
00518 {
00519 long i;
00520 double det,ww,thick;
00521 ivector nodes(nne);
00522 vector x(nne),y(nne),l(3),gp1(intordmm),gp2(intordmm),w(intordmm),t(nne);
00523 matrix n(3,ndofe);
00524
00525 Mt->give_elemnodes (eid,nodes);
00526 Mc->give_thickness (eid,nodes,t);
00527 Mt->give_node_coord2d (x,y,eid);
00528
00529 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00530
00531 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00532
00533 fillm (0.0,lm);
00534
00535 for (i=0;i<intordmm;i++){
00536 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00537 ww=w[i];
00538
00539 bf_matrix (n,x,y,l);
00540
00541 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
00542
00543 ww*=det*thick;
00544
00545 nnj (lm.a,n.a,ww,3,ndofe);
00546 }
00547 }
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560 void cctelem::ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00561 {
00562 long i,ii,ipp;
00563 vector gp1,gp2,w,eps(tncomp),areacoord(3);
00564 matrix gm(tncomp,ndofe);
00565
00566
00567 for (ii=0;ii<nb;ii++){
00568
00569 allocv (intordsm[ii][ii],gp1);
00570 allocv (intordsm[ii][ii],gp2);
00571 allocv (intordsm[ii][ii],w);
00572
00573 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00574
00575 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00576
00577 for (i=0;i<intordsm[ii][ii];i++){
00578 areacoord[0]=gp1[i];
00579 areacoord[1]=gp2[i];
00580 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00581
00582
00583 geom_matrix (gm,x,y,areacoord);
00584
00585 mxv (gm,r,eps);
00586
00587 Mm->storestrain (lcid,ipp,eps);
00588 ipp++;
00589 }
00590 destrv (w); destrv (gp1); destrv (gp2);
00591 }
00592 }
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602 void cctelem::res_ip_strains (long lcid,long eid)
00603 {
00604 vector aux,x(nne),y(nne),r(ndofe);
00605 ivector nodes(nne);
00606 matrix tmat;
00607
00608 Mt->give_node_coord2d (x,y,eid);
00609 Mt->give_elemnodes (eid,nodes);
00610 eldispl (lcid,eid,r.a);
00611
00612
00613 long transf = Mt->locsystems (nodes);
00614 if (transf>0){
00615 allocv (ndofe,aux);
00616 allocm (ndofe,ndofe,tmat);
00617 transf_matrix (nodes,tmat);
00618
00619 lgvectortransf (aux,r,tmat);
00620 copyv (aux,r);
00621 destrv (aux);
00622 destrm (tmat);
00623 }
00624
00625 ip_strains (lcid,eid,0,0,x,y,r);
00626
00627 }
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638 void cctelem::nod_strains_ip (long lcid,long eid,long ri,long ci)
00639 {
00640 long i,j,ipp;
00641 ivector ipnum(nne),nod(nne);
00642 vector eps(tncomp);
00643
00644
00645
00646 ipp=Mt->elements[eid].ipp[ri][ci];
00647 nodip_planelt (ipp,intordsm[0][0],ipnum);
00648
00649
00650 Mt->give_elemnodes (eid,nod);
00651
00652 for (i=0;i<nne;i++){
00653
00654 Mm->givestrain (lcid,ipnum[i],eps);
00655
00656
00657 j=nod[i];
00658 Mt->nodes[j].storestrain (lcid,0,eps);
00659 }
00660
00661 }
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671 void cctelem::res_ip_stresses (long lcid,long eid)
00672 {
00673 long ri,ci;
00674 ri=0;
00675 ci=0;
00676
00677 compute_nlstress (lcid,eid,ri,ci);
00678 }
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689 void cctelem::nod_stresses_ip (long lcid,long eid,long ri,long ci)
00690 {
00691 long i,j,ipp;
00692 ivector ipnum(nne),nod(nne);
00693 vector sig(tncomp);
00694
00695
00696
00697 ipp=Mt->elements[eid].ipp[ri][ci];
00698 nodip_planelt (ipp,intordsm[0][0],ipnum);
00699
00700
00701 Mt->give_elemnodes (eid,nod);
00702
00703 for (i=0;i<nne;i++){
00704
00705 Mm->givestress (lcid,ipnum[i],sig);
00706
00707
00708 j=nod[i];
00709 Mt->nodes[j].storestress (lcid,0,sig);
00710 }
00711
00712 }
00713
00714
00715 void cctelem::strains (long lcid,long eid,long ri,long ci)
00716 {
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774 }
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784 void cctelem::nodecoord (vector &xi,vector &eta)
00785 {
00786 xi[0] = 0.0; eta[0] = 0.0;
00787 xi[1] = 1.0; eta[1] = 0.0;
00788 xi[2] = 0.0; eta[2] = 1.0;
00789 }
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800 void cctelem::appval (double xi,double eta,long fi,long nc,vector &eps,double **val)
00801 {
00802 long i,j,k;
00803 vector nodval(nne),areacoord(3);
00804
00805 areacoord[0]=xi; areacoord[1]=eta; areacoord[2]=1.0-areacoord[0]-areacoord[1];
00806 k=0;
00807 for (i=fi;i<fi+nc;i++){
00808 for (j=0;j<nne;j++){
00809 nodval[j]=val[j][i];
00810 }
00811 eps[k]=approx (areacoord,nodval);
00812 k++;
00813 }
00814 }
00815
00816
00817 void cctelem::stresses (long lcid,long eid,long ri,long ci)
00818 {
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881 }
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897 void cctelem::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00898 {
00899 integratedquant iq;
00900 iq=locstress;
00901
00902
00903 compute_nlstress (lcid,eid,ri,ci);
00904
00905 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00906 }
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921 void cctelem::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00922 {
00923 integratedquant iq;
00924 iq=nonlocstress;
00925
00926
00927 compute_nonloc_nlstress (lcid,eid,ri,ci);
00928
00929 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00930 }
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945 void cctelem::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00946 {
00947 integratedquant iq;
00948 iq=stressincr;
00949
00950
00951 compute_nlstressincr (lcid,eid,ri,ci);
00952
00953 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00954 }
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968 void cctelem::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y)
00969 {
00970 integratedquant iq;
00971 iq=eigstress;
00972
00973
00974 if (Mp->eigstrcomp)
00975 compute_eigstress (lcid,eid,ri,ci);
00976
00977
00978 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y);
00979 }
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992 void cctelem::res_internal_forces (long lcid,long eid,vector &ifor)
00993 {
00994 long transf;
00995 ivector nodes (nne);
00996 vector v(ndofe),x(nne),y(nne);
00997
00998 Mt->give_node_coord2d (x,y,eid);
00999 internal_forces (lcid,eid,0,0,ifor,x,y);
01000
01001
01002 Mt->give_elemnodes (eid,nodes);
01003 transf = Mt->locsystems (nodes);
01004 if (transf>0){
01005 matrix tmat (ndofe,ndofe);
01006 transf_matrix (nodes,tmat);
01007 glvectortransf (ifor,v,tmat);
01008 copyv (v,ifor);
01009 }
01010 }
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023 void cctelem::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
01024 {
01025 long transf;
01026 ivector nodes (nne);
01027 vector v(ndofe),x(nne),y(nne);
01028
01029 Mt->give_node_coord2d (x,y,eid);
01030 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y);
01031
01032
01033 Mt->give_elemnodes (eid,nodes);
01034 transf = Mt->locsystems (nodes);
01035 if (transf>0){
01036 matrix tmat (ndofe,ndofe);
01037 transf_matrix (nodes,tmat);
01038 glvectortransf (ifor,v,tmat);
01039 copyv (v,ifor);
01040 }
01041 }
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054 void cctelem::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01055 {
01056 long transf;
01057 ivector nodes (nne);
01058 vector v(ndofe),x(nne),y(nne);
01059
01060 Mt->give_node_coord2d (x,y,eid);
01061 incr_internal_forces (lcid,eid,0,0,ifor,x,y);
01062
01063
01064 Mt->give_elemnodes (eid,nodes);
01065 transf = Mt->locsystems (nodes);
01066 if (transf>0){
01067 matrix tmat (ndofe,ndofe);
01068 transf_matrix (nodes,tmat);
01069 glvectortransf (ifor,v,tmat);
01070 copyv (v,ifor);
01071 }
01072 }
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085 void cctelem::res_eigstrain_forces (long lcid,long eid,vector &nfor)
01086 {
01087 long transf;
01088 ivector nodes (nne);
01089 vector v(ndofe),x(nne),y(nne);
01090
01091 Mt->give_node_coord2d (x,y,eid);
01092 eigstrain_forces (lcid,eid,0,0,nfor,x,y);
01093
01094
01095 Mt->give_elemnodes (eid,nodes);
01096 transf = Mt->locsystems (nodes);
01097 if (transf>0){
01098 matrix tmat (ndofe,ndofe);
01099 transf_matrix (nodes,tmat);
01100 glvectortransf (nfor,v,tmat);
01101 copyv (v,nfor);
01102 }
01103 }
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116 void cctelem::compute_nlstress (long lcid,long eid,long ri,long ci)
01117 {
01118 long i,ii,ipp;
01119 double thick;
01120 ivector nodes(nne);
01121 vector sig(5),gp1,gp2,w,l(3),t(nne);
01122
01123 Mt->give_elemnodes (eid,nodes);
01124 Mc->give_thickness (eid,nodes,t);
01125
01126 for (ii=0;ii<nb;ii++){
01127 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01128 allocv (intordsm[ii][ii],w);
01129 allocv (intordsm[ii][ii],gp1);
01130 allocv (intordsm[ii][ii],gp2);
01131 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[0][0]);
01132
01133 for (i=0;i<intordsm[ii][ii];i++){
01134
01135 if (Mp->strcomp==1){
01136 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
01137 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
01138
01139 Mm->computenlstresses (ipp);
01140
01141 Mm->givestress (lcid,ipp,sig);
01142 sig[0]*=thick*thick*thick;
01143 sig[1]*=thick*thick*thick;
01144 sig[2]*=thick*thick*thick;
01145 sig[3]*=thick*5.0/6.0;
01146 sig[4]*=thick*5.0/6.0;
01147 Mm->storestress (lcid,ipp,sig);
01148
01149 ipp++;
01150 }
01151 }
01152 destrv (gp1);
01153 destrv (gp2);
01154 destrv (w);
01155 }
01156 }
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170 void cctelem::local_values(long lcid,long eid,long ri,long ci)
01171 {
01172 long i,ii,ipp;
01173
01174 for (ii=0;ii<nb;ii++){
01175 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01176
01177 for (i=0;i<intordsm[ii][ii];i++){
01178
01179 if (Mp->strcomp==1)
01180 Mm->computenlstresses (ipp);
01181 ipp++;
01182 }
01183 }
01184 }
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197 void cctelem::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01198 {
01199 long i,ii,ipp;
01200
01201 for (ii=0;ii<nb;ii++){
01202 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01203
01204 for (i=0;i<intordsm[ii][ii];i++){
01205
01206 if (Mp->strcomp==1)
01207 Mm->compnonloc_nlstresses (ipp);
01208 ipp++;
01209 }
01210 }
01211 }
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224 void cctelem::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01225 {
01226 long i,ii,ipp;
01227
01228 for (ii=0;ii<nb;ii++){
01229 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01230
01231 for (i=0;i<intordsm[ii][ii];i++){
01232
01233 if (Mp->strcomp==1)
01234 Mm->computenlstressesincr (ipp);
01235 ipp++;
01236 }
01237 }
01238 }
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251 void cctelem::compute_eigstress(long lcid,long eid,long ri,long ci)
01252 {
01253 long i,ii,ipp;
01254 vector eigstr(tncomp),sig(tncomp);
01255 matrix d(tncomp,tncomp);
01256
01257 for (ii=0;ii<nb;ii++){
01258 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01259
01260 for (i=0;i<intordsm[ii][ii];i++){
01261 Mm->giveeigstrain (ipp,eigstr);
01262
01263 Mm->matstiff (d,ipp);
01264 mxv (d,eigstr,sig);
01265 Mm->storeeigstress (ipp,sig);
01266 ipp++;
01267 }
01268 }
01269 }
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286 void cctelem::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y)
01287 {
01288 long i,ii,ipp;
01289 double jac;
01290 ivector nodes(nne);
01291 vector w,gp1,gp2,areacoord(3);
01292 vector ipv,contr(ndofe);
01293 matrix gm;
01294
01295 Mt->give_elemnodes (eid,nodes);
01296 fillv (0.0,nv);
01297
01298 for (ii=0;ii<nb;ii++){
01299 allocv (intordsm[ii][ii],gp1);
01300 allocv (intordsm[ii][ii],gp2);
01301 allocv (intordsm[ii][ii],w);
01302 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01303 ipp=Mt->elements[eid].ipp[ii][ii];
01304
01305 for (i=0;i<intordsm[ii][ii];i++){
01306 areacoord[0]=gp1[i];
01307 areacoord[1]=gp2[i];
01308 areacoord[2]=1.0-areacoord[0]-areacoord[2];
01309
01310 jac_2d (jac,x,y,areacoord[0],areacoord[1]);
01311 allocv (ncomp[ii],ipv);
01312
01313
01314 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
01315
01316
01317 allocm (ncomp[ii],ndofe,gm);
01318 geom_matrix_block (gm,ii,x,y,areacoord);
01319
01320 mtxv (gm,ipv,contr);
01321 destrm (gm);
01322 cmulv (jac*w[i],contr);
01323
01324
01325 addv(contr,nv,nv);
01326 ipp++;
01327 destrv (ipv);
01328 }
01329 destrv (w); destrv (gp2); destrv (gp1);
01330 }
01331 }
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357 void cctelem::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
01358 {
01359 long i, j, k, ipp;
01360 long ii, jj, nv = nodval.n;
01361 double xi, eta, ipval;
01362 vector w, gp1, gp2, anv(nne);
01363 long nstra, nstre, ncompstr, ncompeqother;
01364 long idstra, idstre, idoth, idic;
01365 inictype ict;
01366
01367 nstra = idstra = nstre = idstre = idoth = idic = 0;
01368
01369 ict = ictn[0];
01370 for (i=0; i<nne; i++)
01371 {
01372 if (ictn[i] != ict)
01373 {
01374 print_err("Incompatible types of initial conditions on element %ld\n"
01375 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
01376 abort();
01377 }
01378 }
01379 for (j = 0; j < nv; j++)
01380 {
01381 for(i = 0; i < nne; i++)
01382 anv[i] = nodval[i][j];
01383 for (ii = 0; ii < nb; ii++)
01384 {
01385 for (jj = 0; jj < nb; jj++)
01386 {
01387 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01388 if (intordsm[ii][jj] == 0)
01389 continue;
01390 allocv (intordsm[ii][jj],gp1);
01391 allocv (intordsm[ii][jj],gp2);
01392 allocv (intordsm[ii][jj],w);
01393 gauss_points_tr (gp1.a, gp2.a, w.a, intordsm[ii][jj]);
01394 for (k = 0; k < intordsm[ii][jj]; k++)
01395 {
01396 xi=gp1[k];
01397 eta=gp2[k];
01398
01399 ipval = approx_nat (xi,eta,anv);
01400 ncompstr = Mm->ip[ipp].ncompstr;
01401 ncompeqother = Mm->ip[ipp].ncompeqother;
01402 if ((ictn[0] & inistrain) && (j < ncompstr))
01403 {
01404 Mm->ip[ipp].strain[idstra] += ipval;
01405 ipp++;
01406 continue;
01407 }
01408 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01409 {
01410 Mm->ip[ipp].stress[idstre] += ipval;
01411 ipp++;
01412 continue;
01413 }
01414 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
01415 {
01416 Mm->ip[ipp].eqother[idoth] += ipval;
01417 ipp++;
01418 continue;
01419 }
01420 if ((ictn[0] & inicond) && (j < nv))
01421 {
01422 if (Mm->ic[ipp] == NULL)
01423 {
01424 Mm->ic[ipp] = new double[nv-j];
01425 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
01426 }
01427 Mm->ic[ipp][idic] += ipval;
01428 ipp++;
01429 continue;
01430 }
01431 ipp++;
01432 }
01433 destrv(gp1); destrv (gp2); destrv (w);
01434 }
01435 }
01436 ipp=Mt->elements[eid].ipp[ri][ci];
01437 ncompstr = Mm->ip[ipp].ncompstr;
01438 ncompeqother = Mm->ip[ipp].ncompeqother;
01439 if ((ictn[0] & inistrain) && (j < ncompstr))
01440 {
01441 nstra++;
01442 idstra++;
01443 continue;
01444 }
01445 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01446 {
01447 nstre++;
01448 idstre++;
01449 continue;
01450 }
01451 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
01452 {
01453 idoth++;
01454 continue;
01455 }
01456 if ((ictn[0] & inicond) && (j < nv))
01457 {
01458 idic++;
01459 continue;
01460 }
01461 }
01462 }
01463
01464
01465