00001 #include "soilplatetr.h"
00002 #include "dkt.h"
00003 #include "global.h"
00004 #include "globmat.h"
00005 #include "genfile.h"
00006 #include "node.h"
00007 #include "element.h"
00008 #include <stdlib.h>
00009 #include <math.h>
00010
00011 soilplatetr::soilplatetr (void)
00012 {
00013 long i,j;
00014
00015 nne=3; ndofe=9; tncomp=5; napfun=3; ned=3; nned=2;
00016 intordmm=3; intordb=2; ssst=plates;
00017 nb=1;
00018
00019 ncomp = new long [nb];
00020 ncomp[0]=3;
00021
00022 cncomp = new long [nb];
00023 cncomp[0]=0;
00024
00025 nip = new long* [nb];
00026 intordsm = new long* [nb];
00027 for (i=0;i<nb;i++){
00028 nip[i] = new long [nb];
00029 intordsm[i] = new long [nb];
00030 }
00031 nip[0][0]=7;
00032
00033 tnip=0;
00034 for (i=0;i<nb;i++){
00035 for (j=0;j<nb;j++){
00036 tnip+=nip[i][j];
00037 }
00038 }
00039
00040 intordsm[0][0]=7;
00041 }
00042
00043 soilplatetr::~soilplatetr (void)
00044 {
00045 long i;
00046
00047 for (i=0;i<nb;i++){
00048 delete [] nip[i];
00049 delete [] intordsm[i];
00050 }
00051 delete [] nip;
00052 delete [] intordsm;
00053
00054 delete [] cncomp;
00055 delete [] ncomp;
00056 }
00057
00058
00059 double soilplatetr::approx_nat (double xi,double eta,vector &nodval)
00060 {
00061 double f;
00062 vector areacoord(3);
00063
00064
00065 areacoord[0]=xi;
00066 areacoord[1]=eta;
00067 areacoord[2]=1.0-areacoord[0]-areacoord[1];
00068
00069 scprd (areacoord,nodval,f);
00070
00071 return f;
00072 }
00073
00074
00075
00076
00077
00078 void soilplatetr::dbmat (matrix &db,double c1, double c2)
00079 {
00080
00081 db[0][0] = c1; db[0][1] = 0.0; db[0][2] = 0.0;
00082 db[1][0] = 0.0; db[1][1] = c2; db[1][2] = 0.0;
00083 db[2][0] = 0.0; db[2][1] = 0.0; db[2][2] = c2;
00084 }
00085
00086
00087 void soilplatetr::geom_matrix (matrix &gm,matrix &cn,vector &ax,vector &ay,vector &dl, vector &l)
00088 {
00089 int i,i1,i2,i4;
00090 double lll;
00091
00092
00093
00094
00095
00096 lll=l(0)*l(1)*l(2)/2.;
00097 for (i=0;i<nne;i++){
00098 i1=i+1;
00099 if (i1>=nne) i1=i1-nne;
00100 i2=i1+1;
00101 if (i2>=nne) i2=i2-nne;
00102 i4=i*3;
00103 gm[0][i4] = l[i];
00104 gm[0][i4+1]=(l[i]*l[i]*l[i1]+lll)*dl[i]*cn[i][0]-(l[i]*l[i]*l[i2]+lll)*dl[i2]*cn[i2][0];
00105 gm[0][i4+2]=(l[i]*l[i]*l[i1]+lll)*dl[i]*cn[i][1]-(l[i]*l[i]*l[i2]+lll)*dl[i2]*cn[i2][1];
00106 gm[1][i4] = ay[i];
00107
00108
00109
00110
00111
00112
00113
00114
00115 gm[2][i4] = ax[i];
00116
00117
00118
00119
00120
00121
00122
00123
00124 }
00125 }
00126
00127 void soilplatetr::transf_matrix (ivector &nodes,matrix &tmat)
00128 {
00129 long i,n,m;
00130 n=nodes.n;
00131 m=tmat.m;
00132 for (i=0;i<m;i++){
00133 tmat[i][i]=1.0;
00134 }
00135
00136 for (i=0;i<n;i++){
00137 if (Mt->nodes[nodes[i]].transf>0){
00138 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];
00139 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];
00140 }
00141 }
00142 }
00143 void soilplatetr::stiffness_matrix (long eid,long ri, long ci, matrix &sm,vector &x, vector &y)
00144 {
00145 long i,ipp;
00146 double det,jac,thick;
00147 ivector nodes(nne);
00148 vector w,gp1,gp2,c1(nne),c2(nne),t(nne),l(3),dl(3),ax(3),ay(3);
00149 matrix gm,cn(3,2),cc(3,3);
00150
00151 Mt->give_elemnodes (eid,nodes);
00152 Mc->give_thickness (eid,nodes,t);
00153
00154
00155
00156
00157
00158 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]));
00159 dl[0] = sqrt((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
00160 dl[1] = sqrt((x[2]-x[1])*(x[2]-x[1])+(y[2]-y[1])*(y[2]-y[1]));
00161 dl[2] = sqrt((x[0]-x[2])*(x[0]-x[2])+(y[0]-y[2])*(y[0]-y[2]));
00162
00163 cn[0][0]= (y[1]-y[0])/dl[0];
00164 cn[0][1]=-(x[1]-x[0])/dl[0];
00165 cn[1][1]=-(x[2]-x[1])/dl[1];
00166 cn[1][0]= (y[2]-y[1])/dl[1];
00167 cn[2][1]=-(x[0]-x[2])/dl[2];
00168 cn[2][0]= (y[0]-y[2])/dl[2];
00169
00170 ax[0]=(x[2]-x[1])/det; ax[1]=(x[0]-x[2])/det; ax[2]=(x[1]-x[0])/det;
00171
00172 ay[0]=(y[1]-y[2])/det; ay[1]=(y[2]-y[0])/det; ay[2]=(y[0]-y[1])/det;
00173
00174 fillm (0.0,sm);
00175
00176
00177 allocv (intordsm[0][0],w); allocv (intordsm[0][0],gp1); allocv (intordsm[0][0],gp2);
00178 allocm (ncomp[0],ndofe,gm); allocm (ncomp[0],ncomp[0],cc);
00179 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[0][0]);
00180 ipp=Mt->elements[eid].ipp[ri][ci];
00181 for (i=0;i<intordsm[0][0];i++){
00182 thick = approx_nat (gp1[i],gp2[i],t);
00183 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00184 jac=det*w[i];
00185 geom_matrix (gm,cn,ax,ay,dl,l);
00186 Mm->matstiff (cc,ipp);
00187
00188
00189
00190 bdbjac (sm,gm,cc,gm,jac*thick);
00191 ipp++;
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 destrm (cc); destrm (gm); destrv (gp1);destrv (gp2); destrv (w);
00207 }
00208
00209 void soilplatetr::res_stiffness_matrix (long eid,matrix &sm)
00210 {
00211 vector x(nne),y(nne);
00212 Mt->give_node_coord2d (x,y,eid);
00213 stiffness_matrix (eid,0,0,sm,x,y);
00214 }
00215
00216 void soilplatetr::res_mainip_strains (long lcid,long eid)
00217 {
00218 vector aux,x(nne),y(nne),r(ndofe);
00219 ivector nodes(nne);
00220 matrix tmat;
00221
00222 Mt->give_node_coord2d (x,y,eid);
00223 Mt->give_elemnodes (eid,nodes);
00224 eldispl (lcid,eid,r.a);
00225
00226
00227 long transf = Mt->locsystems (nodes);
00228 if (transf>0){
00229 allocv (ndofe,aux);
00230 allocm (ndofe,ndofe,tmat);
00231 transf_matrix (nodes,tmat);
00232
00233 lgvectortransf (aux,r,tmat);
00234 copyv (aux,r);
00235 destrv (aux);
00236 destrm (tmat);
00237 }
00238
00239 mainip_strains (lcid,eid,0,0,x,y,r);
00240 }
00241
00242
00243 void soilplatetr::mainip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00244 {
00245 long i,ipp;
00246 double det;
00247 vector ax(3),ay(3),dl(3),l(3),gp1,gp2,w,eps;
00248 matrix cn(3,2),gm;
00249
00250
00251
00252 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
00253 dl[0] = sqrt((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
00254 dl[1] = sqrt((x[2]-x[1])*(x[2]-x[1])+(y[2]-y[1])*(y[2]-y[1]));
00255 dl[2] = sqrt((x[0]-x[2])*(x[0]-x[2])+(y[0]-y[2])*(y[0]-y[2]));
00256
00257 cn[0][0]= (y[1]-y[0])/dl[0];
00258 cn[0][1]=-(x[1]-x[0])/dl[0];
00259 cn[1][1]=-(x[2]-x[1])/dl[1];
00260 cn[1][0]= (y[2]-y[1])/dl[1];
00261 cn[2][1]=-(x[0]-x[2])/dl[2];
00262 cn[2][0]= (y[0]-y[2])/dl[2];
00263
00264 ax[0]=(x[2]-x[1])/2./det; ax[1]=(x[0]-x[2])/2./det; ax[2]=(x[1]-x[0])/2./det;
00265
00266 ay[0]=(y[1]-y[2])/2./det; ay[1]=(y[2]-y[0])/2./det; ay[2]=(y[0]-y[1])/2./det;
00267
00268
00269
00270 ipp=Mt->elements[eid].ipp[ri][ci];
00271 allocv (intordsm[0][0],w); allocv (intordsm[0][0],gp1); allocv (intordsm[0][0],gp2);
00272 allocm (ncomp[0],ndofe,gm);
00273 allocv (ncomp[0],eps);
00274 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[0][0]);
00275
00276 for (i=0;i<intordsm[0][0];i++){
00277 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00278 geom_matrix (gm,cn,ax,ay,dl,l);
00279 mxv (gm,r,eps);
00280 Mm->storestrain (lcid,ipp,cncomp[0],ncomp[0],eps);
00281 ipp++;
00282
00283
00284 }
00285 destrm (gm); destrv (eps); destrv (gp1); destrv (gp2); destrv (w);
00286
00287 }
00288
00289 void soilplatetr::res_allip_stresses (long lcid,long eid)
00290 {
00291 allip_stresses (lcid, eid, 0, 0);
00292 }
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 void soilplatetr::allip_stresses (long lcid,long eid,long ri,long ci)
00304 {
00305 long ipp,i;
00306 ivector nodes(nne);
00307 vector eps, sig;
00308 matrix cc;
00309
00310 Mt->give_elemnodes (eid,nodes);
00311
00312
00313
00314
00315 allocv (ncomp[0],eps); allocv (ncomp[0],sig);
00316 allocm (ncomp[0],ncomp[0],cc);
00317
00318 ipp=Mt->elements[eid].ipp[ri][ci];
00319
00320
00321 for (i=0;i<intordsm[0][0];i++){
00322 Mm->matstiff (cc,ipp);
00323 Mm->givestrain (lcid,ipp,cncomp[0],ncomp[0],eps);
00324 mxv (cc,eps,sig);
00325 Mm->storestress (lcid,ipp,cncomp[0],sig);
00326
00327
00328
00329 ipp++;
00330 }
00331 destrv (eps); destrv (sig); destrm (cc);
00332
00333
00334
00335
00336 }
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 void soilplatetr::compute_nlstress (long lcid,long eid,long ri,long ci)
00349 {
00350 long i,ii,jj,ipp;
00351
00352 for (ii=0;ii<nb;ii++){
00353 for (jj=0;jj<nb;jj++){
00354 if (intordsm[ii][jj]==0) continue;
00355
00356 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00357
00358 for (i=0;i<intordsm[ii][jj];i++)
00359 {
00360
00361 if (Mp->strcomp==1)
00362 Mm->computenlstresses (ipp);
00363
00364 ipp++;
00365 }
00366 }
00367 }
00368 }
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380 void soilplatetr::compute_nlstressincr (long lcid,long eid,long ri,long ci)
00381 {
00382 long i,ii,jj,ipp;
00383
00384 for (ii=0;ii<nb;ii++){
00385 for (jj=0;jj<nb;jj++){
00386 if (intordsm[ii][jj]==0) continue;
00387
00388 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00389
00390 for (i=0;i<intordsm[ii][jj];i++){
00391
00392 if (Mp->strcomp==1)
00393 Mm->computenlstressesincr (ipp);
00394 ipp++;
00395 }
00396 }
00397 }
00398 }
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 void soilplatetr::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00413 {
00414 integratedquant iq;
00415 iq=locstress;
00416
00417
00418 compute_nlstress (lcid,eid,ri,ci);
00419
00420
00421 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00422 }
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 void soilplatetr::res_internal_forces (long lcid,long eid,vector &ifor)
00434 {
00435 long transf;
00436 ivector nodes (nne);
00437 vector v(ndofe),x(nne),y(nne);
00438
00439 Mt->give_node_coord2d (x,y,eid);
00440
00441 internal_forces (lcid,eid,0,0,ifor,x,y);
00442
00443
00444
00445 Mt->give_elemnodes (eid,nodes);
00446 transf = Mt->locsystems (nodes);
00447 if (transf>0){
00448 matrix tmat (ndofe,ndofe);
00449 transf_matrix (nodes,tmat);
00450
00451 glvectortransf (ifor,v,tmat);
00452 copyv (v,ifor);
00453 }
00454 }
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468 void soilplatetr::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00469 {
00470 integratedquant iq;
00471 iq=stressincr;
00472
00473
00474 compute_nlstressincr (lcid,eid,ri,ci);
00475
00476
00477 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00478 }
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489 void soilplatetr::res_incr_internal_forces (long lcid,long eid,vector &ifor)
00490 {
00491 long transf;
00492 ivector nodes (nne);
00493 vector v(ndofe),x(nne),y(nne);
00494
00495 Mt->give_node_coord2d (x,y,eid);
00496
00497 incr_internal_forces (lcid,eid,0,0,ifor,x,y);
00498
00499
00500
00501 Mt->give_elemnodes (eid,nodes);
00502 transf = Mt->locsystems (nodes);
00503 if (transf>0){
00504 matrix tmat (ndofe,ndofe);
00505 transf_matrix (nodes,tmat);
00506
00507 glvectortransf (ifor,v,tmat);
00508 copyv (v,ifor);
00509 }
00510 }
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526 void soilplatetr::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y)
00527 {
00528 long i,ii,ipp;
00529 double thick,det;
00530 ivector nodes(nne);
00531 vector w,gp1,gp2,t(nne),ipv,contr(ndofe),l(3),dl(3),ax(3),ay(3);
00532 matrix gm,cn(3,2);
00533
00534 Mt->give_elemnodes (eid,nodes);
00535 Mc->give_thickness (eid,nodes,t);
00536
00537 fillv (0.0,nv);
00538
00539
00540 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]));
00541 dl[0] = sqrt((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
00542 dl[1] = sqrt((x[2]-x[1])*(x[2]-x[1])+(y[2]-y[1])*(y[2]-y[1]));
00543 dl[2] = sqrt((x[0]-x[2])*(x[0]-x[2])+(y[0]-y[2])*(y[0]-y[2]));
00544
00545 cn[0][0]= (y[1]-y[0])/dl[0];
00546 cn[0][1]=-(x[1]-x[0])/dl[0];
00547 cn[1][1]=-(x[2]-x[1])/dl[1];
00548 cn[1][0]= (y[2]-y[1])/dl[1];
00549 cn[2][1]=-(x[0]-x[2])/dl[2];
00550 cn[2][0]= (y[0]-y[2])/dl[2];
00551
00552 ax[0]=(x[2]-x[1])/det; ax[1]=(x[0]-x[2])/det; ax[2]=(x[1]-x[0])/det;
00553
00554 ay[0]=(y[1]-y[2])/det; ay[1]=(y[2]-y[0])/det; ay[2]=(y[0]-y[1])/det;
00555
00556
00557
00558 for (ii=0;ii<nb;ii++){
00559 allocv (intordsm[ii][ii],gp1);
00560 allocv (intordsm[ii][ii],gp2);
00561 allocv (intordsm[ii][ii],w);
00562 allocm (ncomp[ii],ndofe,gm);
00563 allocv (ncomp[ii],ipv);
00564
00565 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00566 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00567
00568 for (i=0;i<intordsm[ii][ii];i++){
00569 thick = approx_nat (gp1[i],gp2[i],t);
00570
00571 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
00572
00573
00574 l[0]=gp1[i];
00575 l[1]=gp2[i];
00576 l[2]=1.0-l[0]-l[1];
00577 geom_matrix (gm,cn,ax,ay,dl,l);
00578
00579
00580 mtxv (gm,ipv,contr);
00581
00582 cmulv (det*w[i]*thick,contr);
00583
00584
00585 addv(contr,nv,nv);
00586
00587 ipp++;
00588 }
00589 destrv (ipv); destrm (gm); destrv (w); destrv (gp1); destrv (gp2);
00590 }
00591 }