00001 #include "soilplateq.h"
00002 #include "q4plate.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 soilplateq::soilplateq (void)
00012 {
00013 long i,j;
00014
00015 nne=4; ndofe=12; tncomp=5; napfun=2; ned=4; nned=2;
00016 intordmm=3; intordb=2; ssst=plates;
00017 nb=2;
00018
00019 ncomp = new long [nb];
00020 ncomp[0]=3;
00021 ncomp[1]=2;
00022
00023 cncomp = new long [nb];
00024 cncomp[0]=0;
00025 cncomp[1]=3;
00026
00027 nip = new long* [nb];
00028 intordsm = new long* [nb];
00029 for (i=0;i<nb;i++){
00030 nip[i] = new long [nb];
00031 intordsm[i] = new long [nb];
00032 }
00033 nip[0][0]=9; nip[0][1]=0;
00034 nip[1][0]=0; nip[1][1]=4;
00035
00036 tnip=0;
00037 for (i=0;i<nb;i++){
00038 for (j=0;j<nb;j++){
00039 tnip+=nip[i][j];
00040 }
00041 }
00042
00043 intordsm[0][0]=3; intordsm[0][1]=0;
00044 intordsm[1][0]=0; intordsm[1][1]=2;
00045 }
00046
00047 soilplateq::~soilplateq (void)
00048 {
00049 long i;
00050
00051 for (i=0;i<nb;i++){
00052 delete [] nip[i];
00053 delete [] intordsm[i];
00054 }
00055 delete [] nip;
00056 delete [] intordsm;
00057
00058 delete [] cncomp;
00059 delete [] ncomp;
00060 }
00061
00062
00063
00064
00065
00066
00067
00068 double soilplateq::approx (double xi,double eta,vector &nodval)
00069 {
00070 double f;
00071 vector bf(nne);
00072
00073 bf_lin_4_2d (bf.a,xi,eta);
00074
00075 scprd (bf,nodval,f);
00076
00077 return f;
00078 }
00079
00080
00081
00082
00083
00084 void soilplateq::atd_matrix (matrix &atd,vector &x,vector &y)
00085 {
00086 double sx,sy;
00087 matrix p(16,16),p1(16,16),ct(4,2);
00088 vector dl(4);
00089 long i,j,i1,i2,i3,i4,ii;
00090
00091 dl[0]=sqrt((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
00092 dl[1]=sqrt((x[2]-x[1])*(x[2]-x[1])+(y[2]-y[1])*(y[2]-y[1]));
00093 dl[2]=sqrt((x[3]-x[2])*(x[3]-x[2])+(y[3]-y[2])*(y[3]-y[2]));
00094 dl[3]=sqrt((x[0]-x[3])*(x[0]-x[3])+(y[0]-y[3])*(y[0]-y[3]));
00095
00096 ct[0][0]=(x[1]-x[0])/dl[0];
00097 ct[0][1]=(y[1]-y[0])/dl[0];
00098 ct[1][0]=(x[2]-x[1])/dl[1];
00099 ct[1][1]=(y[2]-y[1])/dl[1];
00100 ct[2][0]=(x[3]-x[2])/dl[2];
00101 ct[2][1]=(y[3]-y[2])/dl[2];
00102 ct[3][0]=(x[0]-x[3])/dl[3];
00103 ct[3][1]=(y[0]-y[3])/dl[3];
00104 fillm (0.0,p);
00105 fillm (0.0,atd);
00106
00107 for (i=0;i<nne;i++){
00108 i1=3*i;
00109 i2=i1+1;
00110 i3=i1+2;
00111 i4=12+i;
00112 p[i1][0]=1.;
00113 p[i1][1]=x[i];
00114 p[i1][2]=y[i];
00115 p[i1][3]=x[i]*y[i];
00116 p[i1][4]=x[i]*x[i];
00117 p[i1][5]=y[i]*y[i];
00118 p[i1][6]=x[i]*x[i]*y[i];
00119 p[i1][7]=y[i]*y[i]*x[i];
00120 p[i1][8]=x[i]*x[i]*x[i];
00121 p[i1][9]=y[i]*y[i]*y[i];
00122 p[i1][10]=x[i]*x[i]*x[i]*y[i];
00123 p[i1][11]=y[i]*y[i]*y[i]*x[i];
00124 p[i2][5]=y[i]*2.;
00125 p[i2][6]=x[i]*x[i];
00126 p[i2][7]=y[i]*x[i]*2.;
00127 p[i2][9]=y[i]*y[i]*3.;
00128 p[i2][10]=x[i]*x[i]*x[i];
00129 p[i2][11]=y[i]*y[i]*x[i]*3.;
00130 p[i2][14]=1.;
00131 p[i2][15]=x[i];
00132 p[i3][4]=-x[i]*2.;
00133 p[i3][6]=-x[i]*y[i]*2.;
00134 p[i3][7]=-y[i]*y[i];
00135 p[i3][8]=-x[i]*x[i]*3.;
00136 p[i3][10]=-x[i]*x[i]*y[i]*3.;
00137 p[i3][11]=-y[i]*y[i]*y[i];
00138 p[i3][12]=-1.;
00139 p[i3][13]=-y[i];
00140
00141 ii=i+1;
00142
00143 if (ii > 3) ii=ii-4;
00144
00145 sy=(y[i]+y[ii])/2.;
00146 sx=(x[i]+x[ii])/2.;
00147 p[i4][1]= ct[i][0];
00148 p[i4][2]= ct[i][1];
00149 p[i4][3]= sy*ct[i][0]+sx*ct[i][1];
00150 p[i4][12]=-ct[i][0];
00151 p[i4][13]=-sy*ct[i][0];
00152 p[i4][14]=-ct[i][1];
00153 p[i4][15]=-sx*ct[i][1];
00154 }
00155
00156 double *ee;
00157 ee = new double [16*16];
00158 for (i=0;i<16;i++){
00159 for (j=0;j<16;j++){
00160 ee[i*16+j]=0.0;
00161 }
00162 ee[i*16+i]=1.0;
00163 }
00164
00165 gemp (p.a,p1.a,ee,16,16,Mp->zero,1);
00166 delete [] ee;
00167
00168 for (i=0;i<p1.n;i++){
00169 atd[i][0] =p1[i][0] -p1[i][12]/dl[0]+p1[i][15]/dl[3];
00170 atd[i][1] =p1[i][1] -p1[i][12]*ct[0][1]/2.-p1[i][15]*ct[3][1]/2.;
00171 atd[i][2] =p1[i][2] +p1[i][12]*ct[0][0]/2.+p1[i][15]*ct[3][0]/2.;
00172 atd[i][3] =p1[i][3] -p1[i][13]/dl[1]+p1[i][12]/dl[0];
00173 atd[i][4] =p1[i][4] -p1[i][13]*ct[1][1]/2.-p1[i][12]*ct[0][1]/2.;
00174 atd[i][5] =p1[i][5] +p1[i][13]*ct[1][0]/2.+p1[i][12]*ct[0][0]/2.;
00175 atd[i][6] =p1[i][6] -p1[i][14]/dl[2]+p1[i][13]/dl[1];
00176 atd[i][7] =p1[i][7] -p1[i][14]*ct[2][1]/2.-p1[i][13]*ct[1][1]/2.;
00177 atd[i][8] =p1[i][8] +p1[i][14]*ct[2][0]/2.+p1[i][13]*ct[1][0]/2.;
00178 atd[i][9] =p1[i][9] -p1[i][15]/dl[3]+p1[i][14]/dl[2];
00179 atd[i][10]=p1[i][10] -p1[i][15]*ct[3][1]/2.-p1[i][14]*ct[2][1]/2.;
00180 atd[i][11]=p1[i][11] +p1[i][15]*ct[3][0]/2.+p1[i][14]*ct[2][0]/2.;
00181 }
00182
00183 }
00184
00185
00186
00187
00188
00189 void soilplateq::dbmat (matrix &db,double c1, double c2, long ri, long ci)
00190 {
00191
00192 if (ri==0 && ci==0){
00193 db[0][0] = c1; db[0][1] = 0.0; db[0][2] = 0.0;
00194 db[1][0] = 0.0; db[1][1] = c2; db[1][2] = 0.0;
00195 db[2][0] = 0.0; db[2][1] = 0.0; db[2][2] = c2;
00196 }
00197 if (ri==1 && ci==1){
00198 db[0][0] = c2*5.0/6.0; db[0][1] = 0.0;
00199 db[1][0] = 0.0; db[1][1] = c2*5.0/6.0;
00200 }
00201 }
00202 void soilplateq::transfmat (ivector &nodes,matrix &tmat)
00203 {
00204 long i,n,m;
00205 n=nodes.n;
00206 m=tmat.m;
00207 for (i=0;i<m;i++){
00208 tmat[i][i]=1.0;
00209 }
00210
00211 for (i=0;i<n;i++){
00212 if (Mt->nodes[nodes[i]].transf>0){
00213 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];
00214 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];
00215 }
00216 }
00217 }
00218
00219 void soilplateq::appval (double xi,double eta,long fi,long nc,vector &eps,double **val)
00220 {
00221 long i,j,k;
00222 vector nodval;
00223
00224 k=0;
00225 allocv (nne,nodval);
00226 for (i=fi;i<fi+nc;i++){
00227 for (j=0;j<nne;j++){
00228 nodval[j]=val[j][i];
00229 }
00230 eps[k]=approx (xi,eta,nodval);
00231 k++;
00232 }
00233
00234 destrv (nodval);
00235 }
00236
00237 void soilplateq::nodecoord (vector &xi,vector &eta)
00238 {
00239 xi[0] = 1.0; eta[0] = 1.0;
00240 xi[1] =-1.0; eta[1] = 1.0;
00241 xi[2] =-1.0; eta[2] =-1.0;
00242 xi[3] = 1.0; eta[3] =-1.0;
00243 }
00244
00245
00246 void soilplateq::geom_matrix (matrix &gm,matrix &atd,vector &x,vector &y, vector &l)
00247 {
00248 int i;
00249 double sx,sy;
00250
00251
00252 sx=(x[0]*(1+l[0])*(1+l[1])+x[1]*(1-l[0])*(1+l[1])+x[2]*(1-l[0])*(1-l[1])+x[3]*(1+l[0])*(1-l[1]))/4.;
00253 sy=(y[0]*(1+l[0])*(1+l[1])+y[1]*(1-l[0])*(1+l[1])+y[2]*(1-l[0])*(1-l[1])+y[3]*(1+l[0])*(1-l[1]))/4.;
00254
00255 for (i=0;i<gm.n;i++){
00256 gm[0][i]=atd[0][i]+sx*atd[1][i]+sy*atd[2][i]+sx*sy*atd[3][i]+sx*sx*atd[4][i]+sy*sy*atd[5][i]+
00257 sx*sx*sy*atd[6][i]+sy*sy*sx*atd[7][i]+sx*sx*sx*atd[8][i]+sy*sy*sy*atd[9][i]+sx*sx*sx*sy*atd[10][i]+sx*sy*sy*sy*atd[11][i];
00258 gm[1][i]=atd[1][i]+sy*atd[3][i]+2.*sx*atd[4][i]+2.*sx*sy*atd[6][i]+sy*sy*atd[7][i]+3.*sx*sx*atd[8][i]+3.*sx*sx*sy*atd[10][i]+sy*sy*sy*atd[11][i];
00259
00260 gm[2][i]=atd[2][i]+sx*atd[3][i]+2.*sy*atd[5][i]+sx*sx*atd[6][i]+2.*sy*sx*atd[7][i]+3.*sy*sy*atd[9][i]+sx*sx*sx*atd[10][i]+3.*sx*sy*sy*atd[11][i];
00261
00262 }
00263 }
00264 void soilplateq::geom_matrix_shear (matrix &gm,matrix &atd,vector &x,vector &y,vector &l)
00265 {
00266 long i;
00267 double sx,sy;
00268
00269
00270 sx=(x[0]*(1+l[0])*(1+l[1])+x[1]*(1-l[0])*(1+l[1])+x[2]*(1-l[0])*(1-l[1])+x[3]*(1+l[0])*(1-l[1]))/4.;
00271 sy=(y[0]*(1+l[0])*(1+l[1])+y[1]*(1-l[0])*(1+l[1])+y[2]*(1-l[0])*(1-l[1])+y[3]*(1+l[0])*(1-l[1]))/4.;
00272
00273 for (i=0;i<atd.n;i++){
00274 gm[0][i]=atd[1][i]+sy*atd[3][i]-atd[12][i]-sy*atd[13][i];
00275 gm[1][i]=atd[2][i]+sx*atd[3][i]-atd[14][i]-sx*atd[15][i];
00276 }
00277 }
00278
00279 void soilplateq::res_stiffness_matrix (long eid,matrix &sm)
00280 {
00281 vector x(nne),y(nne);
00282 Mt->give_node_coord2d (x,y,eid);
00283 stiffness_matrix (eid,0,0,sm,x,y);
00284 }
00285
00286
00287 void soilplateq::stiffness_matrix (long eid,long ri, long ci, matrix &sm,vector &x, vector &y)
00288 {
00289 long ii,jj,i,j,ipp;
00290 double jac,a,a0,a1;
00291 ivector nodes(nne);
00292 vector w,gp,l(2);
00293 matrix gm,cc,atd(16,12),d(3,3);
00294
00295 Mt->give_elemnodes (eid,nodes);
00296
00297
00298 a=(x[0]+x[1]+x[2]+x[3])/4.;
00299 x[0]=x[0]-a;
00300 x[1]=x[1]-a;
00301 x[2]=x[2]-a;
00302 x[3]=x[3]-a;
00303 a=(y[0]+y[1]+y[2]+y[3])/4.;
00304 y[0]=y[0]-a;
00305 y[1]=y[1]-a;
00306 y[2]=y[2]-a;
00307 y[3]=y[3]-a;
00308
00309
00310
00311 a = (x[2]-x[0])*(y[3]-y[1])-(y[2]-y[0])*(x[3]-x[1]);
00312 a0=-(x[3]-x[2])*(y[1]-y[0])+(y[3]-y[2])*(x[1]-x[0]);
00313 a1=-(x[2]-x[1])*(y[3]-y[0])+(y[2]-y[1])*(x[3]-x[0]);
00314
00315 atd_matrix (atd,x,y);
00316 fillm (0.0,sm);
00317
00318 ipp=Mt->elements[eid].ipp[ri][ci];
00319 Mm->matstiff (d,ipp);
00320
00321 for (ii=0;ii<nb;ii++){
00322 for (jj=0;jj<nb;jj++){
00323 if (intordsm[ii][jj]==0) continue;
00324 allocv (intordsm[ii][jj],w); allocv (intordsm[ii][jj],gp);
00325 allocm (ncomp[jj],ndofe,gm);
00326 allocm (ncomp[ii],ncomp[jj],cc);
00327
00328 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00329 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00330
00331 for (i=0;i<intordsm[ii][jj];i++){
00332 l[0]=gp[i];
00333 for (j=0;j<intordsm[ii][jj];j++){
00334 l[1]=gp[j];
00335 jac=(a+a0*l[0]+a1*l[1])/8.*w[i]*w[j];
00336
00337 if(ii==0) geom_matrix (gm,atd,x,y,l);
00338 else if(ii==1) geom_matrix_shear (gm,atd,x,y,l);
00339 dbmat (cc, d[0][0],d[1][1],ii,jj);
00340 bdbjac (sm,gm,cc,gm,jac);
00341 ipp++;
00342 }
00343 }
00344 destrm (cc); destrm (gm); destrv (gp); destrv (w);
00345 }
00346
00347 }
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 }
00373
00374
00375 void soilplateq::res_mainip_strains (long lcid,long eid)
00376 {
00377 vector aux,x(nne),y(nne),r(ndofe);
00378 ivector nodes(nne);
00379 matrix tmat;
00380
00381 Mt->give_node_coord2d (x,y,eid);
00382 Mt->give_elemnodes (eid,nodes);
00383 eldispl (lcid,eid,r.a);
00384
00385
00386 long transf = Mt->locsystems (nodes);
00387 if (transf>0){
00388 allocv (ndofe,aux);
00389 allocm (ndofe,ndofe,tmat);
00390 transfmat (nodes,tmat);
00391
00392 lgvectortransf (aux,r,tmat);
00393 copyv (aux,r);
00394 destrv (aux);
00395 destrm (tmat);
00396 }
00397
00398 mainip_strains (lcid,eid,0,0,x,y,r);
00399 }
00400
00401
00402 void soilplateq::mainip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00403 {
00404 long i,j,ii,ipp;
00405 double a;
00406 vector gp,w,eps,natcoord(2),l(3);
00407 matrix gm,atd(16,12);
00408
00409
00410 a=(x[0]+x[1]+x[2]+x[3])/4.;
00411 x[0]=x[0]-a;
00412 x[1]=x[1]-a;
00413 x[2]=x[2]-a;
00414 x[3]=x[3]-a;
00415 a=(y[0]+y[1]+y[2]+y[3])/4.;
00416 y[0]=y[0]-a;
00417 y[1]=y[1]-a;
00418 y[2]=y[2]-a;
00419 y[3]=y[3]-a;
00420 atd_matrix (atd,x,y);
00421
00422
00423 for (ii=0;ii<nb;ii++){
00424 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
00425 allocv (intordsm[ii][ii],gp); allocv (intordsm[ii][ii],w);
00426 allocv (ncomp[ii],eps);
00427 allocm (ncomp[ii],ndofe,gm);
00428 gauss_points (gp.a,w.a,intordsm[ii][ii]);
00429
00430 for (i=0;i<intordsm[ii][ii];i++){
00431 l[0]=gp[i];
00432 for (j=0;j<intordsm[ii][ii];j++){
00433 l[1]=gp[j];
00434 if(ii==0) geom_matrix (gm,atd,x,y,l);
00435 else if(ii==1) geom_matrix_shear (gm,atd,x,y,l);
00436 mxv (gm,r,eps);
00437 Mm->storestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00438 ipp++;
00439
00440 }
00441 }
00442 destrm (gm); destrv (eps); destrv (gp);
00443 }
00444
00445 }
00446
00447 void soilplateq::elem_strains (double **stra,long lcid,long eid,long ri,long ci)
00448 {
00449 long i,j,ii,ipp;
00450 double *lsm,*lhs,*rhs;
00451 vector nxi(nne),neta(nne),gp,w,eps,natcoord(2),l(3);
00452 ivector cn,nodes(nne);
00453
00454 lsm = new double [9];
00455 nullv (lsm,9);
00456
00457 nodecoord (nxi,neta);
00458 Mt->give_elemnodes (eid,nodes);
00459
00460
00461
00462 for (ii=0;ii<nb;ii++){
00463 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
00464 allocv (ncomp[ii],eps);
00465 allocv (intordsm[ii][ii],gp); allocv (intordsm[ii][ii],w);
00466 lhs = new double [ncomp[ii]*3];
00467 rhs = new double [ncomp[ii]*3];
00468 gauss_points (gp.a,w.a,intordsm[ii][ii]);
00469 nullv (rhs,ncomp[ii]*3);
00470
00471 for (i=0;i<intordsm[ii][ii];i++){
00472 l[0]=gp[i];
00473 for (j=0;j<intordsm[ii][ii];j++){
00474 l[1]=gp[j];
00475 Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00476 natcoord[0]=l[0]; natcoord[1]=l[1];
00477 matassem_lsm (lsm,natcoord);
00478 rhsassem_lsm (rhs,natcoord,eps);
00479
00480
00481
00482
00483 ipp++;
00484 }
00485 }
00486 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
00487 nodal_values (stra,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii]);
00488 delete [] lhs; delete [] rhs;
00489 destrv (eps); destrv (w); destrv (gp);
00490 }
00491 delete [] lsm;
00492 }
00493
00494 void soilplateq::nod_strains (long lcid,long eid,long ri,long ci)
00495 {
00496 long i,j,ii,ipp;
00497 double *lsm,*lhs,*rhs;
00498 vector nxi(nne),neta(nne),gp,w,eps,natcoord(2),l(3);
00499 ivector cn,nodes(nne);
00500
00501 lsm = new double [9];
00502
00503 nodecoord (nxi,neta);
00504 Mt->give_elemnodes (eid,nodes);
00505
00506
00507
00508
00509 for (ii=0;ii<nb;ii++){
00510 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
00511 allocv (ncomp[ii],eps);
00512 allocv (intordsm[ii][ii],gp); allocv (intordsm[ii][ii],w);
00513 lhs = new double [ncomp[ii]*3];
00514 rhs = new double [ncomp[ii]*3];
00515 gauss_points (gp.a,w.a,intordsm[ii][ii]);
00516 nullv (lsm,9);
00517 nullv (rhs,ncomp[ii]*3);
00518
00519 for (i=0;i<intordsm[ii][ii];i++){
00520 l[0]=gp[i];
00521 for (j=0;j<intordsm[ii][ii];j++){
00522 l[1]=gp[j];
00523 Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00524 natcoord[0]=l[0]; natcoord[1]=l[1];
00525 matassem_lsm (lsm,natcoord);
00526 rhsassem_lsm (rhs,natcoord,eps);
00527
00528
00529
00530
00531 ipp++;
00532 }
00533 }
00534 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
00535 Mt->strain_nodal_values (nodes,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii],lcid);
00536 delete [] lhs; delete [] rhs;
00537 destrv (eps); destrv (w); destrv (gp);
00538 }
00539
00540 delete [] lsm;
00541 }
00542
00543 void soilplateq::appstrain (long lcid,long eid,double xi,double eta,long fi,long ncomp,vector &eps)
00544 {
00545 long i,j,k;
00546 ivector nodes;
00547 vector l(3),nodval;
00548 if (ncomp != eps.n){
00549 fprintf (stderr,"\n\n wrong interval of indices in function strain (%s, line %d).\n",__FILE__,__LINE__);
00550 abort ();
00551 }
00552
00553
00554 allocv (nne,nodes);
00555 allocv (nne,nodval);
00556 Mt->give_elemnodes (eid,nodes);
00557 k=0;
00558 for (i=fi;i<fi+ncomp;i++){
00559 for (j=0;j<nne;j++){
00560 nodval[j]=Mt->nodes[nodes[j]].strain[lcid*tncomp+i];
00561 }
00562 eps[k]=approx (xi,eta,nodval);
00563 k++;
00564 }
00565
00566 destrv (nodes); destrv (nodval);
00567 }
00568
00569 void soilplateq::allip_strains (double **stra,long lcid,long eid,long ri,long ci)
00570 {
00571 long i,j,ii,jj,ipp;
00572 vector eps(tncomp),gp,w,areacoord(3);
00573
00574 for (ii=0;ii<nb;ii++){
00575 for (jj=0;jj<nb;jj++){
00576 if (intordsm[ii][jj]==0) continue;
00577
00578 allocv (intordsm[ii][jj],gp); allocv (intordsm[ii][jj],w);
00579 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00580 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00581
00582 for (i=0;i<intordsm[ii][jj];i++){
00583 for (j=0;j<intordsm[ii][jj];j++){
00584 if (Mp->strainaver==0)
00585 appval (gp[i],gp[j],0,tncomp,eps,stra);
00586 if (Mp->strainaver==1)
00587 appstrain (lcid,eid,gp[i],gp[j],0,tncomp,eps);
00588 Mm->storestrain (lcid,ipp,eps);
00589 ipp++;
00590 }
00591 }
00592 destrv (w); destrv (gp);
00593 }
00594 }
00595 }
00596
00597 void soilplateq::strains (long lcid,long eid,long ri,long ci)
00598 {
00599 long i,naep,ncp,sid;
00600 double **stra;
00601 vector coord,eps;
00602
00603 if (Mp->strainaver==0){
00604 stra = new double* [nne];
00605 for (i=0;i<nne;i++){
00606 stra[i] = new double [tncomp];
00607 }
00608 elem_strains (stra,lcid,eid,ri,ci);
00609 }
00610
00611 switch (Mm->stra.tape[eid]){
00612 case nowhere:{
00613 break;
00614 }
00615 case intpts:{
00616 allip_strains (stra,lcid,eid,ri,ci);
00617 break;
00618 }
00619 case enodes:{
00620 break;
00621 }
00622 case userdefined:{
00623
00624 naep = Mm->stra.give_naep (eid);
00625 ncp = Mm->stra.give_ncomp (eid);
00626 sid = Mm->stra.give_sid (eid);
00627 allocv (ncp,eps);
00628 allocv (2,coord);
00629 for (i=0;i<naep;i++){
00630 Mm->stra.give_aepcoord (sid,i,coord);
00631 if (Mp->strainaver==0)
00632 appval (coord[0],coord[1],0,ncp,eps,stra);
00633 if (Mp->strainaver==1)
00634 appstrain (lcid,eid,coord[0],coord[1],0,ncp,eps);
00635 Mm->stra.storevalues(lcid,eid,i,eps);
00636 }
00637 destrv (eps);
00638 destrv (coord);
00639 break;
00640 }
00641 default:{
00642 fprintf (stderr,"\n\n unknown strain point is required in function dstelem::strains (%s, line %d).\n",__FILE__,__LINE__);
00643 }
00644 }
00645 if (Mp->strainaver==0){
00646 for (i=0;i<nne;i++){
00647 delete [] stra[i];
00648 }
00649 delete [] stra;
00650 }
00651
00652 }
00653
00654
00655 void soilplateq::res_allip_stresses (long lcid,long eid)
00656 {
00657 allip_stresses (lcid, eid, 0, 0);
00658 }
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669 void soilplateq::allip_stresses (long lcid,long eid,long ri,long ci)
00670 {
00671 long ipp,i,j;
00672 ivector nodes(nne);
00673 vector eps, sig;
00674 matrix cc;
00675
00676 Mt->give_elemnodes (eid,nodes);
00677
00678
00679
00680
00681 allocv (ncomp[0],eps); allocv (ncomp[0],sig);
00682 allocm (ncomp[0],ncomp[0],cc);
00683
00684 ipp=Mt->elements[eid].ipp[ri][ci];
00685
00686
00687 for (i=0;i<intordsm[0][0];i++){
00688 for (j=0;j<intordsm[0][0];j++){
00689 Mm->matstiff (cc,ipp);
00690 Mm->givestrain (lcid,ipp,cncomp[0],ncomp[0],eps);
00691 mxv (cc,eps,sig);
00692 Mm->storestress (lcid,ipp,cncomp[0],sig);
00693
00694
00695
00696 ipp++;
00697 }
00698 }
00699 destrv (eps); destrv (sig); destrm (cc);
00700
00701
00702
00703
00704 }
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722 void soilplateq::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y,vector &z)
00723 {
00724 integratedquant iq;
00725 iq=locstress;
00726
00727
00728 compute_nlstress (lcid,eid,ri,ci);
00729
00730
00731 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
00732 }
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751 void soilplateq::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y,vector &z)
00752 {
00753 integratedquant iq;
00754 iq=nonlocstress;
00755
00756
00757 compute_nonloc_nlstress (lcid,eid,ri,ci);
00758
00759
00760 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
00761 }
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776 void soilplateq::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y,vector &z)
00777 {
00778 integratedquant iq;
00779 iq=stressincr;
00780
00781
00782 compute_nlstressincr (lcid,eid,ri,ci);
00783
00784
00785 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y,z);
00786 }
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 void soilplateq::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y,vector &z)
00806 {
00807 integratedquant iq;
00808 iq=eigstress;
00809
00810
00811 if (Mp->eigstrcomp)
00812 compute_eigstress (lcid,eid,ri,ci);
00813
00814
00815 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y,z);
00816 }
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829 void soilplateq::res_internal_forces (long lcid,long eid,vector &ifor)
00830 {
00831
00832 ivector nodes (nne);
00833 vector v(ndofe),x(nne),y(nne),z(nne);
00834
00835 Mt->give_node_coord3d (x,y,z,eid);
00836
00837 internal_forces (lcid,eid,0,0,ifor,x,y,z);
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 void soilplateq::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
00866 {
00867
00868 ivector nodes (nne);
00869 vector v(ndofe),x(nne),y(nne),z(nne);
00870
00871 Mt->give_node_coord3d (x,y,z,eid);
00872
00873 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y,z);
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888 }
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901 void soilplateq::res_incr_internal_forces (long lcid,long eid,vector &ifor)
00902 {
00903
00904 ivector nodes (nne);
00905 vector v(ndofe),x(nne),y(nne),z(nne);
00906
00907 Mt->give_node_coord3d (x,y,z,eid);
00908
00909 incr_internal_forces (lcid,eid,0,0,ifor,x,y,z);
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924 }
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938 void soilplateq::res_eigstrain_forces (long lcid,long eid,vector &nfor)
00939 {
00940
00941 ivector nodes (nne);
00942 vector v(ndofe),x(nne),y(nne),z(nne);
00943
00944 Mt->give_node_coord3d (x,y,z,eid);
00945
00946 eigstrain_forces (lcid,eid,0,0,nfor,x,y,z);
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961 }
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974 void soilplateq::compute_nlstress (long lcid,long eid,long ri,long ci)
00975 {
00976 long ipp,ii,jj,i,j;
00977
00978 for (ii=0;ii<nb;ii++){
00979 for (jj=0;jj<nb;jj++){
00980 if (intordsm[ii][jj]==0) continue;
00981
00982 ipp=Mt->elements[eid].ipp[ii+ri][jj+ci];
00983
00984 for (i=0;i<intordsm[ii][jj];i++){
00985 for (j=0;j<intordsm[ii][jj];j++){
00986
00987 if (Mp->strcomp==1)
00988 {
00989 Mm->computenlstresses (ipp);
00990 ipp++;
00991 }
00992 }
00993 }
00994 }
00995 }
00996 }
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010 void soilplateq::local_values (long lcid,long eid,long ri,long ci)
01011 {
01012 long ipp,ii,i,j;
01013 vector eps(3),sig(3);
01014 matrix cc;
01015
01016 for (ii=0;ii<nb;ii++){
01017 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
01018 allocm (ncomp[ii],ncomp[ii],cc);
01019 for (i=0;i<intordsm[0][0];i++){
01020 for (j=0;j<intordsm[0][0];j++){
01021
01022 if (Mp->strcomp==1)
01023 {
01024 Mm->computenlstresses (ipp);
01025 Mm->givestress (lcid,ipp,eps);
01026 Mm->matstiff (cc,ipp);
01027 mtxv (cc,eps,sig);
01028 Mm->storestress (lcid,ipp,sig);
01029 }
01030 ipp++;
01031 }
01032 }
01033 destrm (cc);
01034 }
01035 }
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048 void soilplateq::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01049 {
01050 long ipp,ii,i,j;
01051 vector eps(3),sig(3);
01052 matrix cc;
01053
01054 for (ii=0;ii<nb;ii++){
01055 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
01056 allocm (ncomp[ii],ncomp[ii],cc);
01057 for (i=0;i<intordsm[0][0];i++){
01058 for (j=0;j<intordsm[0][0];j++){
01059
01060 if (Mp->strcomp==1)
01061 {
01062 Mm->compnonloc_nlstresses (ipp);
01063 Mm->givestress (lcid,ipp,eps);
01064 Mm->matstiff (cc,ipp);
01065 mtxv (cc,eps,sig);
01066 Mm->storestress (lcid,ipp,sig);
01067 }
01068 ipp++;
01069 }
01070 }
01071 destrm (cc);
01072 }
01073 }
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086 void soilplateq::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01087 {
01088 long ipp,ii,jj,i,j;
01089
01090 for (ii=0;ii<nb;ii++){
01091 for (jj=0;jj<nb;jj++){
01092 if (intordsm[ii][jj]==0) continue;
01093
01094 ipp=Mt->elements[eid].ipp[ii+ri][jj+ci];
01095
01096 for (i=0;i<intordsm[ii][jj];i++){
01097 for (j=0;j<intordsm[ii][jj];j++){
01098
01099 if (Mp->strcomp==1)
01100 {
01101 Mm->computenlstressesincr (ipp);
01102 ipp++;
01103 }
01104 }
01105 }
01106 }
01107 }
01108 }
01109
01110
01111
01112 void soilplateq::compute_eigstress(long lcid,long eid,long ri,long ci)
01113 {
01114 long ipp,ii,i,j;
01115 vector eigstr(tncomp),sig(tncomp);
01116 matrix d(tncomp,tncomp);
01117
01118 for (ii=0;ii<nb;ii++){
01119 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
01120 for (i=0;i<intordsm[0][0];i++){
01121 for (j=0;j<intordsm[0][0];j++){
01122 Mm->giveeigstrain (ipp,eigstr);
01123
01124 Mm->matstiff (d,ipp);
01125 mxv (d,eigstr,sig);
01126 Mm->storeeigstress (ipp,sig);
01127 ipp++;
01128 }
01129 }
01130 }
01131 }
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148 void soilplateq::elem_integration (integratedquant iq, long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y, vector&z)
01149 {
01150 long ipp,ii,i,j;
01151 double a,a0,a1,jac;
01152 ivector nodes(nne);
01153 vector f(ndofe),ipv(3);
01154 vector w,gp,l(2);
01155 matrix gm,cc,tran(3,3),atd(16,12);
01156
01157 Mt->give_elemnodes (eid,nodes);
01158 Mt->give_node_coord3d (x,y,z,eid);
01159
01160
01161 a=(x[0]+x[1]+x[2]+x[3])/4.;
01162 x[0]=x[0]-a;
01163 x[1]=x[1]-a;
01164 x[2]=x[2]-a;
01165 x[3]=x[3]-a;
01166 a=(y[0]+y[1]+y[2]+y[3])/4.;
01167 y[0]=y[0]-a;
01168 y[1]=y[1]-a;
01169 y[2]=y[2]-a;
01170 y[3]=y[3]-a;
01171
01172 a = (x[2]-x[0])*(y[3]-y[1])-(y[2]-y[0])*(x[3]-x[1]);
01173 a0=-(x[3]-x[2])*(y[1]-y[0])+(y[3]-y[2])*(x[1]-x[0]);
01174 a1=-(x[2]-x[1])*(y[3]-y[0])+(y[2]-y[1])*(x[3]-x[0]);
01175 atd_matrix (atd,x,y);
01176 fillv (0.0,nv);
01177
01178 for (ii=0;ii<nb;ii++){
01179 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
01180 allocv (intordsm[ii][ii],gp); allocv (intordsm[ii][ii],w);
01181 allocm (ncomp[ii],ndofe,gm);
01182 gauss_points (gp.a,w.a,intordsm[ii][ii]);
01183
01184 for (i=0;i<intordsm[0][0];i++){
01185 for (j=0;j<intordsm[0][0];j++){
01186 l[0]=gp[i]; l[1]=gp[j];
01187 jac=(a+a0*l[0]+a1*l[1])/8.*w[i]*w[j];
01188
01189 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
01190
01191 geom_matrix (gm,atd,x,y,l);
01192
01193 mtxv (gm,ipv,f);
01194 cmulv (jac,f);
01195
01196 addv(f,nv,nv);
01197 ipp++;
01198 }
01199 }
01200 destrm (gm); destrv (gp); destrv (w);
01201 }
01202 }