00001 #include "dst.h"
00002 #include "global.h"
00003 #include "globmat.h"
00004 #include "genfile.h"
00005 #include "node.h"
00006 #include "element.h"
00007 #include "intpoints.h"
00008 #include <stdlib.h>
00009 #include <math.h>
00010
00011 dstelem::dstelem (void)
00012 {
00013 long i,j;
00014
00015
00016 nne=3;
00017
00018 ndofe=9;
00019
00020 tncomp=5;
00021
00022 napfun=3;
00023
00024 ned=3;
00025
00026 nned=2;
00027 intordmm=3;
00028
00029 intordb=2;
00030
00031 ssst=plates;
00032
00033
00034 nb=3;
00035
00036
00037 ncomp = new long [nb];
00038 ncomp[0]=3;
00039 ncomp[1]=3;
00040 ncomp[2]=2;
00041
00042
00043 cncomp = new long [nb];
00044 cncomp[0]=0;
00045 cncomp[1]=0;
00046 cncomp[2]=3;
00047
00048
00049
00050 nip = new long* [nb];
00051 intordsm = new long* [nb];
00052 for (i=0;i<nb;i++){
00053 nip[i] = new long [nb];
00054 intordsm[i] = new long [nb];
00055 }
00056 nip[0][0]=1; nip[0][1]=0; nip[0][2]=0;
00057 nip[1][0]=0; nip[1][1]=3; nip[1][2]=0;
00058 nip[2][0]=0; nip[2][1]=0; nip[2][2]=1;
00059
00060
00061 tnip=0;
00062 for (i=0;i<nb;i++){
00063 for (j=0;j<nb;j++){
00064 tnip+=nip[i][j];
00065 }
00066 }
00067
00068 intordsm[0][0]=1; intordsm[0][1]=0; intordsm[0][2]=0;
00069 intordsm[1][0]=0; intordsm[1][1]=3; intordsm[1][2]=0;
00070 intordsm[2][0]=0; intordsm[2][1]=0; intordsm[2][2]=1;
00071 }
00072
00073 dstelem::~dstelem (void)
00074 {
00075 long i;
00076
00077 for (i=0;i<nb;i++){
00078 delete [] nip[i];
00079 delete [] intordsm[i];
00080 }
00081 delete [] nip;
00082 delete [] intordsm;
00083
00084 delete [] cncomp;
00085 delete [] ncomp;
00086 }
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 double dstelem::approx (vector &areacoord,vector &nodval)
00098 {
00099 double f;
00100 scprd (areacoord,nodval,f);
00101 return f;
00102 }
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 void dstelem::tran_matrix (matrix &a,matrix &ct,long eid)
00115 {
00116 double det,dj,gg,thick;
00117 long i,ipp;
00118 ivector nodes(nne);
00119 vector x(3),y(3),ax(3),ay(3),dl(3),t(nne);
00120 matrix dbx(3,3),dby(3,3),b(2,3),ba(3,3),dd(3,3),d(tncomp,tncomp);
00121
00122
00123
00124 Mt->give_elemnodes (eid,nodes);
00125 Mt->give_node_coord2d (x,y,eid);
00126 ipp=Mt->elements[eid].ipp[0][0];
00127 Mm->matstiff (d,ipp);
00128 Mc->give_thickness (eid,nodes,t);
00129 thick = (t[0]+t[1]+t[2])/3.0;
00130 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
00131 dl[0] = sqrt((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
00132 dl[1] = sqrt((x[2]-x[1])*(x[2]-x[1])+(y[2]-y[1])*(y[2]-y[1]));
00133 dl[2] = sqrt((x[0]-x[2])*(x[0]-x[2])+(y[0]-y[2])*(y[0]-y[2]));
00134
00135 ct[0][0]=(x[1]-x[0])/dl[0];
00136 ct[0][1]=(y[1]-y[0])/dl[0];
00137 ct[1][0]=(x[2]-x[1])/dl[1];
00138 ct[1][1]=(y[2]-y[1])/dl[1];
00139 ct[2][0]=(x[0]-x[2])/dl[2];
00140 ct[2][1]=(y[0]-y[2])/dl[2];
00141
00142 ax[0]=(x[2]-x[1])/2./det; ax[1]=(x[0]-x[2])/2./det; ax[2]=(x[1]-x[0])/2./det;
00143
00144 ay[0]=(y[1]-y[2])/2./det; ay[1]=(y[2]-y[0])/2./det; ay[2]=(y[0]-y[1])/2./det;
00145
00146 dj=ay[1]*ax[0]+ay[0]*ax[1];
00147 dbx[0][0]= 8.*ct[0][0]*ay[0]*ay[1];
00148 dbx[1][0]= 4.*ct[0][1]*dj;
00149 dbx[2][0]= 4.*ct[0][0]*dj +8.*ct[0][1]*ay[0]*ay[1];
00150 dby[0][0]= 4.*ct[0][0]*dj;
00151 dby[1][0]= 8.*ct[0][1]*ax[0]*ax[1];
00152 dby[2][0]= 8.*ct[0][0]*ax[0]*ax[1] +4.*ct[0][1]*dj;
00153
00154 dj=ay[2]*ax[1]+ay[1]*ax[2];
00155 dbx[0][1]= 8.*ct[1][0]*ay[1]*ay[2];
00156 dbx[1][1]= 4.*ct[1][1]*dj;
00157 dbx[2][1]= 4.*ct[1][0]*dj +8.*ct[1][1]*ay[1]*ay[2];
00158 dby[0][1]= 4.*ct[1][0]*dj;
00159 dby[1][1]= 8.*ct[1][1]*ax[1]*ax[2];
00160 dby[2][1]= 8.*ct[1][0]*ax[1]*ax[2] +4.*ct[1][1]*dj;
00161
00162 dj=ay[0]*ax[2]+ay[2]*ax[0];
00163 dbx[0][2]= 8.*ct[2][0]*ay[2]*ay[0];
00164 dbx[1][2]= 4.*ct[2][1]*dj;
00165 dbx[2][2]= 4.*ct[2][0]*dj +8.*ct[2][1]*ay[2]*ay[0];
00166 dby[0][2]= 4.*ct[2][0]*dj;
00167 dby[1][2]= 8.*ct[2][1]*ax[2]*ax[0];
00168 dby[2][2]= 8.*ct[2][0]*ax[2]*ax[0] +4.*ct[2][1]*dj;
00169
00170 dmatblock (dd, d,0,0,thick);
00171
00172 gg=dd[2][2]/thick/thick*10.;
00173 for (i=0;i<3;i++){
00174 b[0][i]=(dd[0][0]*dbx[0][i]+dd[2][0]*dby[0][i]
00175 +dd[0][1]*dbx[1][i]+dd[2][1]*dby[1][i]
00176 +dd[0][2]*dbx[2][i]+dd[2][2]*dby[2][i])/gg;
00177 b[1][i]=(dd[1][0]*dby[0][i]+dd[2][0]*dbx[0][i]
00178 +dd[1][1]*dby[1][i]+dd[2][1]*dbx[1][i]
00179 +dd[1][2]*dby[2][i]+dd[2][2]*dbx[2][i])/gg;
00180 }
00181
00182
00183 dbx[0][0]=-.5*dl[0];
00184 dbx[1][0]=-dl[1]*(ct[0][0]*ct[1][0]+ct[0][1]*ct[1][1])/6.;
00185 dbx[2][0]=-dl[2]*(ct[0][0]*ct[2][0]+ct[0][1]*ct[2][1])/6.;
00186 dbx[0][1]=-dl[0]*(ct[1][0]*ct[0][0]+ct[1][1]*ct[0][1])/6.;
00187 dbx[1][1]=-.5*dl[1];
00188 dbx[2][1]=-dl[2]*(ct[1][0]*ct[2][0]+ct[1][1]*ct[2][1])/6.;
00189 dbx[0][2]=-dl[0]*(ct[2][0]*ct[0][0]+ct[2][1]*ct[0][1])/6.;
00190 dbx[1][2]=-dl[1]*(ct[2][0]*ct[1][0]+ct[2][1]*ct[1][1])/6.;
00191 dbx[2][2]=-.5*dl[2];
00192
00193 for (i=0;i<3;i++){
00194 dbx[i][0]=dbx[i][0]+dl[i]*(ct[i][0]*b[0][0]+ct[i][1]*b[1][0]);
00195 dbx[i][1]=dbx[i][1]+dl[i]*(ct[i][0]*b[0][1]+ct[i][1]*b[1][1]);
00196 dbx[i][2]=dbx[i][2]+dl[i]*(ct[i][0]*b[0][2]+ct[i][1]*b[1][2]);
00197 }
00198
00199 det= dbx[0][0]*dbx[1][1]*dbx[2][2]+dbx[0][1]*dbx[1][2]*dbx[2][0]+dbx[0][2]*dbx[1][0]*dbx[2][1]
00200 -dbx[0][1]*dbx[1][0]*dbx[2][2]-dbx[0][0]*dbx[1][2]*dbx[2][1]-dbx[0][2]*dbx[1][1]*dbx[2][0];
00201 ba[0][0]=(dbx[1][1]*dbx[2][2]-dbx[1][2]*dbx[2][1])/det;
00202 ba[0][1]=(dbx[0][2]*dbx[2][1]-dbx[0][1]*dbx[2][2])/det;
00203 ba[0][2]=(dbx[0][1]*dbx[1][2]-dbx[0][2]*dbx[1][1])/det;
00204 ba[1][0]=(dbx[1][2]*dbx[2][0]-dbx[1][0]*dbx[2][2])/det;
00205 ba[1][1]=(dbx[0][0]*dbx[2][2]-dbx[0][2]*dbx[2][0])/det;
00206 ba[1][2]=(dbx[0][2]*dbx[1][0]-dbx[0][0]*dbx[1][2])/det;
00207 ba[2][0]=(dbx[1][0]*dbx[2][1]-dbx[1][1]*dbx[2][0])/det;
00208 ba[2][1]=(dbx[0][1]*dbx[2][0]-dbx[0][0]*dbx[2][1])/det;
00209 ba[2][2]=(dbx[0][0]*dbx[1][1]-dbx[0][1]*dbx[1][0])/det;
00210
00211 for (i=0;i<3;i++){
00212 a[i][0]=(-ba[i][0]+ba[i][2]);
00213 a[i][1]=-(ba[i][0]*ct[0][1]*dl[0]+ba[i][2]*ct[2][1]*dl[2])/2.;
00214 a[i][2]= (ba[i][0]*ct[0][0]*dl[0]+ba[i][2]*ct[2][0]*dl[2])/2.;
00215 a[i][3]=( ba[i][0]-ba[i][1]);
00216 a[i][4]=-(ba[i][0]*ct[0][1]*dl[0]+ba[i][1]*ct[1][1]*dl[1])/2.;
00217 a[i][5]= (ba[i][0]*ct[0][0]*dl[0]+ba[i][1]*ct[1][0]*dl[1])/2.;
00218 a[i][6]=( ba[i][1]-ba[i][2]);
00219 a[i][7]=-(ba[i][1]*ct[1][1]*dl[1]+ba[i][2]*ct[2][1]*dl[2])/2.;
00220 a[i][8]= (ba[i][1]*ct[1][0]*dl[1]+ba[i][2]*ct[2][0]*dl[2])/2.;
00221 }
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 }
00232
00233
00234 void dstelem::geom_matrix_bconst (matrix &gm,vector &x,vector &y)
00235 {
00236 double det;
00237 long i,i1,i2;
00238 vector ax(3),ay(3);
00239
00240
00241 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.;
00242
00243 ax[0]=(x[2]-x[1])/2./det; ax[1]=(x[0]-x[2])/2./det; ax[2]=(x[1]-x[0])/2./det;
00244
00245 ay[0]=(y[1]-y[2])/2./det; ay[1]=(y[2]-y[0])/2./det; ay[2]=(y[0]-y[1])/2./det;
00246
00247 for (i=0;i<3;i++){
00248 i2=3*i+2;
00249 i1=i2-1;
00250 gm[0][i2]+=ay[i];
00251 gm[1][i1]+=-ax[i];
00252 gm[2][i2]+=ax[i];
00253 gm[2][i1]+=-ay[i];
00254 }
00255
00256 }
00257
00258 void dstelem::geom_matrix_bending (matrix &gm,matrix &a,matrix &ct,vector &x,vector &y,vector &l)
00259 {
00260 double det;
00261 long i,i2,i3;
00262 vector ax(3),ay(3);
00263 matrix ba(3,3);
00264
00265
00266 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.;
00267
00268 ax[0]=(x[2]-x[1])/2./det; ax[1]=(x[0]-x[2])/2./det; ax[2]=(x[1]-x[0])/2./det;
00269
00270 ay[0]=(y[1]-y[2])/2./det; ay[1]=(y[2]-y[0])/2./det; ay[2]=(y[0]-y[1])/2./det;
00271
00272 for (i=0;i<3;i++){
00273 i2=i+1;
00274 if (i2>2){i2=i2-3;}
00275 i3=i+2;
00276 if (i3>2){i3=i3-3;}
00277 ba[0][i]= 4.*ct[i][0]*( l[i2]*ay[i]+l[i]*ay[i2]+ay[i3]/3.);
00278 ba[1][i]= 4.*ct[i][1]*( l[i2]*ax[i]+l[i]*ax[i2]+ax[i3]/3.);
00279 ba[2][i]= 4.*ct[i][0]*( l[i2]*ax[i]+l[i]*ax[i2]+ax[i3]/3.)+4.*ct[i][1]*( l[i2]*ay[i]+l[i]*ay[i2]+ay[i3]/3.);
00280 }
00281
00282 for (i=0;i<9;i++){
00283 gm[0][i]=ba[0][0]*a[0][i]+ba[0][1]*a[1][i]+ba[0][2]*a[2][i];
00284 gm[1][i]=ba[1][0]*a[0][i]+ba[1][1]*a[1][i]+ba[1][2]*a[2][i];
00285 gm[2][i]=ba[2][0]*a[0][i]+ba[2][1]*a[1][i]+ba[2][2]*a[2][i];
00286
00287 }
00288
00289 }
00290
00291 void dstelem::geom_matrix_shear (matrix &gs,matrix &a,matrix &ct,long eid)
00292 {
00293 double det,dj,gg,thick;
00294 long i,ipp;
00295 ivector nodes(nne);
00296 vector x(3),y(3),ax(3),ay(3),t(nne);
00297 matrix dbx(3,3),dby(3,3),dd(3,3),b(2,3),d(tncomp,tncomp);
00298
00299 Mt->give_elemnodes (eid,nodes);
00300 Mt->give_node_coord2d (x,y,eid);
00301 ipp=Mt->elements[eid].ipp[0][0];
00302 Mm->matstiff (d,ipp);
00303 Mc->give_thickness (eid,nodes,t);
00304 thick = (t[0]+t[1]+t[2])/3.;
00305
00306 det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.;
00307
00308 ax[0]=(x[2]-x[1])/2./det; ax[1]=(x[0]-x[2])/2./det; ax[2]=(x[1]-x[0])/2./det;
00309
00310 ay[0]=(y[1]-y[2])/2./det; ay[1]=(y[2]-y[0])/2./det; ay[2]=(y[0]-y[1])/2./det;
00311 dj=ay[1]*ax[0]+ay[0]*ax[1];
00312 dbx[0][0]= 8.*ct[0][0]*ay[0]*ay[1];
00313 dbx[1][0]= 4.*ct[0][1]*dj;
00314 dbx[2][0]= 4.*ct[0][0]*dj +8.*ct[0][1]*ay[0]*ay[1];
00315 dby[0][0]= 4.*ct[0][0]*dj;
00316 dby[1][0]= 8.*ct[0][1]*ax[0]*ax[1];
00317 dby[2][0]= 8.*ct[0][0]*ax[0]*ax[1] +4.*ct[0][1]*dj;
00318
00319 dj=ay[2]*ax[1]+ay[1]*ax[2];
00320 dbx[0][1]= 8.*ct[1][0]*ay[1]*ay[2];
00321 dbx[1][1]= 4.*ct[1][1]*dj;
00322 dbx[2][1]= 4.*ct[1][0]*dj +8.*ct[1][1]*ay[1]*ay[2];
00323 dby[0][1]= 4.*ct[1][0]*dj;
00324 dby[1][1]= 8.*ct[1][1]*ax[1]*ax[2];
00325 dby[2][1]= 8.*ct[1][0]*ax[1]*ax[2] +4.*ct[1][1]*dj;
00326
00327 dj=ay[0]*ax[2]+ay[2]*ax[0];
00328 dbx[0][2]= 8.*ct[2][0]*ay[2]*ay[0];
00329 dbx[1][2]= 4.*ct[2][1]*dj;
00330 dbx[2][2]= 4.*ct[2][0]*dj +8.*ct[2][1]*ay[2]*ay[0];
00331 dby[0][2]= 4.*ct[2][0]*dj;
00332 dby[1][2]= 8.*ct[2][1]*ax[2]*ax[0];
00333 dby[2][2]= 8.*ct[2][0]*ax[2]*ax[0] +4.*ct[2][1]*dj;
00334
00335 dmatblock (dd, d,0,0,thick);
00336 gg=dd[2][2]/thick/thick*10.;
00337 for (i=0;i<3;i++){
00338 b[0][i]=(dd[0][0]*dbx[0][i]+dd[2][0]*dby[0][i]
00339 +dd[0][1]*dbx[1][i]+dd[2][1]*dby[1][i]
00340 +dd[0][2]*dbx[2][i]+dd[2][2]*dby[2][i])/gg;
00341 b[1][i]=(dd[1][0]*dby[0][i]+dd[2][0]*dbx[0][i]
00342 +dd[1][1]*dby[1][i]+dd[2][1]*dbx[1][i]
00343 +dd[1][2]*dby[2][i]+dd[2][2]*dbx[2][i])/gg;
00344 }
00345 for (i=0;i<9;i++){
00346 gs[0][i]=b[0][0]*a[0][i]+b[0][1]*a[1][i]+b[0][2]*a[2][i];
00347 gs[1][i]=b[1][0]*a[0][i]+b[1][1]*a[1][i]+b[1][2]*a[2][i];
00348 }
00349
00350 }
00351
00352
00353
00354 void dstelem::transf_matrix (ivector &nodes,matrix &tmat)
00355 {
00356 long i,n,m;
00357 n=nodes.n;
00358 m=tmat.m;
00359 for (i=0;i<m;i++){
00360 tmat[i][i]=1.0;
00361 }
00362
00363 for (i=0;i<n;i++){
00364 if (Mt->nodes[nodes[i]].transf>0){
00365 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];
00366 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];
00367 }
00368 }
00369 }
00370
00371
00372
00373
00374
00375
00376 void dstelem::dmatblock (matrix &dd,matrix &d,long ri, long ci, double t)
00377 {
00378 double c;
00379
00380 c=t*t*t;
00381 fillm (0.0,dd);
00382
00383 if ((ri==0 && ci==0)||(ri==1 && ci==1)){
00384 dd[0][0] = c*d[0][0]; dd[0][1] = c*d[0][1]; dd[0][2] = c*d[0][2];
00385 dd[1][0] = c*d[1][0]; dd[1][1] = c*d[1][1]; dd[1][2] = c*d[1][2];
00386 dd[2][0] = c*d[2][0]; dd[2][1] = c*d[2][1]; dd[2][2] = c*d[2][2];
00387 }
00388 if (ri==2 && ci==2){
00389 dd[0][0] = t*d[3][3]*5.0/6.0; dd[0][1] = 0.0;
00390 dd[1][0] = 0.0; dd[1][1] = t*d[4][4]*5.0/6.0;
00391 }
00392 }
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402 void dstelem::stiffness_matrix (long eid,long ri,long ci, matrix &sm,vector &x, vector &y)
00403 {
00404 long i,ii,jj,ipp;
00405 double jac,det,thick;
00406 ivector nodes(nne);
00407 vector w,gp1,gp2,l(3),t(nne);
00408 matrix gmc,dd,a(3,9),ct(3,2),d(tncomp,tncomp);
00409
00410 Mt->give_elemnodes (eid,nodes);
00411 Mc->give_thickness (eid,nodes,t);
00412 tran_matrix (a,ct,eid);
00413
00414
00415 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00416
00417 fillm (0.0,sm);
00418
00419 for (ii=0;ii<nb;ii++){
00420 for (jj=0;jj<nb;jj++){
00421 if (intordsm[ii][jj]==0) continue;
00422
00423 allocv (intordsm[ii][jj],w);
00424 allocv (intordsm[ii][jj],gp1); allocv (intordsm[ii][jj],gp2);
00425 allocm (ncomp[ii],ndofe,gmc);
00426 allocm (ncomp[ii],ncomp[jj],dd);
00427
00428 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
00429
00430 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00431
00432 for (i=0;i<intordsm[ii][jj];i++){
00433 fillm (0.0,gmc);
00434 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00435
00436
00437
00438 if (ii==0){
00439 geom_matrix_bconst (gmc,x,y);
00440 }
00441 if (ii==1){
00442 geom_matrix_bending (gmc,a,ct,x,y,l);
00443 }
00444 if(ii==2){
00445 geom_matrix_shear (gmc,a,ct,eid);
00446 }
00447
00448 Mm->matstiff (d,ipp);
00449 ipp++;
00450 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
00451
00452 dmatblock (dd, d,ii,jj,thick);
00453 jac=det*w[i];
00454 bdbjac (sm,gmc,dd,gmc,jac);
00455
00456
00457
00458 }
00459 destrm (dd); destrm (gmc); destrv (gp1); destrv (gp2); destrv (w);
00460 }
00461 }
00462
00463 }
00464
00465 void dstelem::res_stiffness_matrix (long eid,matrix &sm)
00466 {
00467 long transf;
00468 ivector nodes(nne);
00469 vector x(nne),y(nne);
00470 Mt->give_node_coord2d (x,y,eid);
00471 Mt->give_elemnodes (eid,nodes);
00472
00473 stiffness_matrix (eid,0,0,sm,x,y);
00474
00475
00476 transf = Mt->locsystems (nodes);
00477 if (transf>0){
00478 matrix tmat (ndofe,ndofe);
00479 transf_matrix (nodes,tmat);
00480 glmatrixtransf (sm,tmat);
00481 }
00482 }
00483 void dstelem::nodecoord (vector &xi,vector &eta)
00484 {
00485 xi[0] = 0.0; eta[0] = 0.0;
00486 xi[1] = 1.0; eta[1] = 0.0;
00487 xi[2] = 0.0; eta[2] = 1.0;
00488 }
00489
00490
00491
00492
00493 void dstelem::appval (vector &l, long fi,long nc,vector &eps,double **val)
00494 {
00495 long i,j,k;
00496 vector nodval;
00497
00498 k=0;
00499 allocv (nne,nodval);
00500 for (i=fi;i<fi+nc;i++){
00501 for (j=0;j<nne;j++){
00502 nodval[j]=val[j][i];
00503 }
00504 eps[k] = nodval[0]*l[0]+nodval[1]*l[1]+nodval[2]*l[2];
00505 k++;
00506 }
00507
00508 destrv (nodval);
00509 }
00510
00511
00512
00513
00514
00515
00516 void dstelem::ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00517 {
00518 long ii,i,ipp;
00519 vector eps,w,gp1,gp2,l(3);
00520 matrix gm,a(3,9),ct(3,2);
00521
00522 tran_matrix (a,ct,eid);
00523
00524 for (ii=0;ii<nb;ii++){
00525 ipp=Mt->elements[eid].ipp[ii+ri][ii+ci];
00526 allocm (ncomp[ii],ndofe,gm);
00527 allocv (ncomp[ii],eps);
00528 allocv (intordsm[ii][ii],w);
00529 allocv (intordsm[ii][ii],gp1);
00530 allocv (intordsm[ii][ii],gp2);
00531 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00532
00533 for (i=0;i<intordsm[ii][ii];i++){
00534 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00535
00536 if (ii==0){
00537
00538
00539 geom_matrix_bending (gm,a,ct,x,y,l);
00540 geom_matrix_bconst (gm,x,y);
00541 }
00542 if (ii==1){
00543 geom_matrix_bending (gm,a,ct,x,y,l);
00544 geom_matrix_bconst (gm,x,y);
00545 }
00546 if(ii==2)
00547 geom_matrix_shear (gm,a,ct,eid);
00548
00549 mxv (gm,r,eps);
00550 Mm->storestrain (lcid,ipp,cncomp[ii],eps);
00551 ipp++;
00552 }
00553 destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
00554 destrm (gm);
00555 }
00556 }
00557
00558
00559 void dstelem::res_ip_strains (long lcid,long eid)
00560 {
00561 vector aux,x(nne),y(nne),r(ndofe);
00562 ivector nodes(nne);
00563 matrix tmat;
00564
00565 Mt->give_node_coord2d (x,y,eid);
00566 Mt->give_elemnodes (eid,nodes);
00567 eldispl (lcid,eid,r.a);
00568
00569
00570 long transf = Mt->locsystems (nodes);
00571 if (transf>0){
00572 allocv (ndofe,aux);
00573 allocm (ndofe,ndofe,tmat);
00574 transf_matrix (nodes,tmat);
00575
00576 lgvectortransf (aux,r,tmat);
00577 copyv (aux,r);
00578 destrv (aux);
00579 destrm (tmat);
00580 }
00581
00582 ip_strains (lcid,eid,0,0,x,y,r);
00583
00584 }
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596 void dstelem::res_ip_stresses (long lcid,long eid)
00597 {
00598 compute_nlstress (lcid,eid,0,0);
00599 }
00600
00601 void dstelem::stresses (long lcid,long eid,long ri, long ci)
00602 {
00603
00604 }
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623 void dstelem::compute_nlstress (long lcid,long eid,long ri,long ci)
00624 {
00625 long i,ii,jj,ipp;
00626 double thick;
00627 ivector nodes(nne);
00628 vector sig(5),gp1,gp2,w,l(3),t(nne);
00629
00630 Mt->give_elemnodes (eid,nodes);
00631 Mc->give_thickness (eid,nodes,t);
00632
00633 for (ii=0;ii<nb;ii++){
00634 for (jj=0;jj<nb;jj++){
00635 if (intordsm[ii][jj]==0) continue;
00636
00637 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00638 allocv (intordsm[ii][jj],w);
00639 allocv (intordsm[ii][jj],gp1);
00640 allocv (intordsm[ii][jj],gp2);
00641 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[0][0]);
00642
00643 for (i=0;i<intordsm[ii][jj];i++){
00644
00645 if (Mp->strcomp==1){
00646 l[0]=gp1[i]; l[1]=gp2[i]; l[2]=1.0-l[0]-l[1];
00647 thick = t[0]*l[0]+t[1]*l[1]+t[2]*l[2];
00648
00649 Mm->computenlstresses (ipp);
00650 Mm->givestress (lcid,ipp,sig);
00651 sig[0]*=thick*thick*thick;
00652 sig[1]*=thick*thick*thick;
00653 sig[2]*=thick*thick*thick;
00654 sig[3]*=thick*5.0/6.0;
00655 sig[4]*=thick*5.0/6.0;
00656 Mm->storestress (lcid,ipp,sig);
00657
00658 ipp++;
00659 }
00660 }
00661 destrv (gp1);
00662 destrv (gp2);
00663 destrv (w);
00664 }
00665 }
00666 }
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679 void dstelem::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00680 {
00681 integratedquant iq;
00682 iq=locstress;
00683
00684
00685 compute_nlstress (lcid,eid,ri,ci);
00686
00687 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00688 }
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699 void dstelem::res_internal_forces (long lcid,long eid,vector &ifor)
00700 {
00701 long transf;
00702 ivector nodes (nne);
00703 vector v(ndofe),x(nne),y(nne);
00704
00705 Mt->give_node_coord2d (x,y,eid);
00706 internal_forces (lcid,eid,0,0,ifor,x,y);
00707
00708
00709 Mt->give_elemnodes (eid,nodes);
00710 transf = Mt->locsystems (nodes);
00711 if (transf>0){
00712 matrix tmat (ndofe,ndofe);
00713 transf_matrix (nodes,tmat);
00714 glvectortransf (ifor,v,tmat);
00715 copyv (v,ifor);
00716 }
00717 }
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729 void dstelem::compute_nlstressincr (long lcid,long eid,long ri,long ci)
00730 {
00731 long i,j,ii,jj,ipp;
00732
00733 for (ii=0;ii<nb;ii++){
00734 for (jj=0;jj<nb;jj++){
00735 if (intordsm[ii][jj]==0) continue;
00736
00737 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00738
00739 for (i=0;i<intordsm[ii][jj];i++){
00740 for (j=0;j<intordsm[ii][jj];j++){
00741
00742 if (Mp->strcomp==1)
00743 Mm->computenlstressesincr (ipp);
00744 ipp++;
00745 }
00746 }
00747 }
00748 }
00749 }
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762 void dstelem::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
00763 {
00764 integratedquant iq;
00765 iq=stressincr;
00766
00767
00768 compute_nlstressincr (lcid,eid,ri,ci);
00769
00770
00771 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
00772 }
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783 void dstelem::res_incr_internal_forces (long lcid,long eid,vector &ifor)
00784 {
00785 long transf;
00786 ivector nodes (nne);
00787 vector v(ndofe),x(nne),y(nne);
00788
00789 Mt->give_node_coord2d (x,y,eid);
00790
00791 incr_internal_forces (lcid,eid,0,0,ifor,x,y);
00792
00793
00794
00795 Mt->give_elemnodes (eid,nodes);
00796 transf = Mt->locsystems (nodes);
00797 if (transf>0){
00798 matrix tmat (ndofe,ndofe);
00799 transf_matrix (nodes,tmat);
00800
00801 glvectortransf (ifor,v,tmat);
00802 copyv (v,ifor);
00803 }
00804 }
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819 void dstelem::elem_integration (integratedquant iq,long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y)
00820 {
00821 long i,ii,ipp;
00822 double jac,thick,det;
00823 ivector nodes(nne);
00824 vector w,gp1,gp2,t(nne),ipv,contr(ndofe),l(3);
00825 matrix gm,a(3,9),ct(3,2);
00826
00827 Mt->give_elemnodes (eid,nodes);
00828 Mc->give_thickness (eid,nodes,t);
00829
00830 tran_matrix (a,ct,eid);
00831
00832 det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00833
00834 fillv (0.0,nv);
00835
00836 for (ii=0;ii<nb;ii++){
00837 allocv (intordsm[ii][ii],gp1);
00838 allocv (intordsm[ii][ii],gp2);
00839 allocv (intordsm[ii][ii],w);
00840 allocm (ncomp[ii],ndofe,gm);
00841 allocv (ncomp[ii],ipv);
00842
00843 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00844 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00845
00846 for (i=0;i<intordsm[ii][ii];i++){
00847 l[0]=gp1[i];
00848 l[1]=gp2[i];
00849 l[2]=1.0-l[0]-l[1];
00850
00851 thick = approx (l,t);
00852
00853 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
00854
00855
00856
00857
00858
00859
00860 if(ii==1){
00861 geom_matrix_bending (gm,a,ct,x,y,l);
00862 geom_matrix_bconst (gm,x,y);
00863 }
00864 if(ii==2){
00865 geom_matrix_shear (gm,a,ct,eid);
00866 }
00867
00868
00869 mtxv (gm,ipv,contr);
00870
00871 jac=det*w[i];
00872
00873 cmulv (jac,contr);
00874
00875
00876 addv(contr,nv,nv);
00877
00878 ipp++;
00879 }
00880 destrv (ipv); destrm (gm); destrv (w); destrv (gp1); destrv (gp2);
00881 }
00882 }
00883
00884
00885
00886
00887
00888
00889
00890 void dstelem::nodeforces (long eid,long *le,double *nv,vector &nf)
00891 {
00892 double cn0,sn0,cn1,sn1,cn2,sn2;
00893 ivector nodes(nne);
00894 vector x(nne),y(nne),dl(3);
00895
00896 Mt->give_elemnodes (eid,nodes);
00897 Mt->give_node_coord2d (x,y,eid);
00898
00899 if (le[0]==1){
00900 cn0= (y[1]-y[0])*dl[0]/60.;
00901 sn0=-(x[1]-x[0])*dl[0]/60.;
00902 nf[1]+=-(3.*nv[0]+2.*nv[3])*cn0;
00903 nf[2]+=-(3.*nv[0]+2.*nv[3])*sn0;
00904 nf[4]+= (2.*nv[0]+3.*nv[3])*cn0;
00905 nf[5]+= (2.*nv[0]+3.*nv[3])*sn0;
00906 nf[0]+=((2.-0.1)*nv[0] +(1.+0.1)*nv[3])*dl[0]/6.;
00907 nf[3]+=((1.-0.1)*nv[0] +(2.+0.1)*nv[3])*dl[0]/6.;
00908 }
00909 if (le[1]==1){
00910 cn1= (y[2]-y[1])*dl[1]/60.;
00911 sn1=-(x[2]-x[1])*dl[1]/60.;
00912 nf[4]+=-(3.*nv[3]+2.*nv[6])*cn1;
00913 nf[5]+=-(3.*nv[3]+2.*nv[6])*sn1;
00914 nf[7]+= (2.*nv[3]+3.*nv[6])*cn1;
00915 nf[8]+= (2.*nv[3]+3.*nv[6])*sn1;
00916 nf[3]+=((2.-0.1)*nv[3] +(1.+0.1)*nv[6])*dl[1]/6.;
00917 nf[6]+=((1.-0.1)*nv[3] +(2.+0.1)*nv[6])*dl[1]/6.;
00918 }
00919 if (le[2]==1){
00920 cn2= (y[0]-y[2])*dl[2]/60.;
00921 sn2=-(x[0]-x[2])*dl[2]/60.;
00922 nf[7]+=-(3.*nv[6]+2.*nv[0])*cn2;
00923 nf[8]+=-(3.*nv[6]+2.*nv[0])*sn2;
00924 nf[1]+= (2.*nv[6]+3.*nv[0])*cn2;
00925 nf[2]+= (2.*nv[6]+3.*nv[0])*sn2;
00926 nf[6]+=((2.-0.1)*nv[6] +(1.+0.1)*nv[0])*dl[2]/6.;
00927 nf[0]+=((1.-0.1)*nv[6] +(2.+0.1)*nv[0])*dl[2]/6.;
00928 }
00929 }
00930
00931
00932
00933
00934
00935
00936 void dstelem::areaforces (long eid,double *nv,vector &nf)
00937 {
00938 double pl,cn0,cn1,cn2,sn0,sn1,sn2;
00939 ivector nodes(nne);
00940 vector x(nne),y(nne);
00941
00942 Mt->give_elemnodes (eid,nodes);
00943 Mt->give_node_coord2d (x,y,eid);
00944
00945 pl = fabs(((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.);
00946
00947 cn0= (y[1]-y[0])*pl/360.;
00948 cn1= (y[2]-y[1])*pl/360.;
00949 cn2= (y[0]-y[2])*pl/360.;
00950 sn0=-(x[1]-x[0])*pl/360.;
00951 sn1=-(x[2]-x[1])*pl/360.;
00952 sn2=-(x[0]-x[2])*pl/360.;
00953 nf[1]=-(7.*nv[0]+3.*nv[3]+5.*nv[6])*cn2+(7.*nv[0]+5.*nv[3]+3.*nv[6])*cn0;
00954 nf[2]=-(7.*nv[0]+3.*nv[3]+5.*nv[6])*sn2+(7.*nv[0]+5.*nv[3]+3.*nv[6])*sn0;
00955 nf[4]=-(5.*nv[0]+7.*nv[3]+3.*nv[6])*cn0+(3.*nv[0]+7.*nv[3]+5.*nv[6])*cn1;
00956 nf[5]=-(5.*nv[0]+7.*nv[3]+3.*nv[6])*sn0+(3.*nv[0]+7.*nv[3]+5.*nv[6])*sn1;
00957 nf[7]=-(3.*nv[0]+5.*nv[3]+7.*nv[6])*cn1+(5.*nv[0]+3.*nv[3]+7.*nv[6])*cn2;
00958 nf[8]=-(3.*nv[0]+5.*nv[3]+7.*nv[6])*sn1+(5.*nv[0]+3.*nv[3]+7.*nv[6])*sn2;
00959 nf[0]=(nv[0]/6. +nv[3]/12.+nv[6]/12.)*pl;
00960 nf[3]=(nv[0]/12.+nv[3]/6. +nv[6]/12.)*pl;
00961 nf[6]=(nv[0]/12.+nv[3]/12.+nv[6]/6. )*pl;
00962
00963
00964 }
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990 void dstelem::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
00991 {
00992 long i, j, k, ipp;
00993 long ii, jj, nv = nodval.n;
00994 double xi, eta, ipval;
00995 vector w, gp1, gp2, anv(nne);
00996 long nstra, nstre, ncompstr, ncompeqother;
00997 long idstra, idstre, idoth, idic;
00998 inictype ict;
00999
01000 nstra = idstra = nstre = idstre = idoth = idic = 0;
01001
01002 ict = ictn[0];
01003 for (i=0; i<nne; i++)
01004 {
01005 if (ictn[i] != ict)
01006 {
01007 print_err("Incompatible types of initial conditions on element %ld\n"
01008 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
01009 abort();
01010 }
01011 }
01012 for (j = 0; j < nv; j++)
01013 {
01014 for(i = 0; i < nne; i++)
01015 anv[i] = nodval[i][j];
01016 for (ii = 0; ii < nb; ii++)
01017 {
01018 for (jj = 0; jj < nb; jj++)
01019 {
01020 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01021 if (intordsm[ii][jj] == 0)
01022 continue;
01023 allocv (intordsm[ii][jj],gp1);
01024 allocv (intordsm[ii][jj],gp2);
01025 allocv (intordsm[ii][jj],w);
01026 gauss_points_tr (gp1.a, gp2.a, w.a, intordsm[ii][jj]);
01027 for (k = 0; k < intordsm[ii][jj]; k++)
01028 {
01029 xi=gp1[k];
01030 eta=gp2[k];
01031
01032 ipval = 0.;
01033
01034 ncompstr = Mm->ip[ipp].ncompstr;
01035 ncompeqother = Mm->ip[ipp].ncompeqother;
01036 if ((ictn[0] & inistrain) && (j < ncompstr))
01037 {
01038 Mm->ip[ipp].strain[idstra] += ipval;
01039 ipp++;
01040 continue;
01041 }
01042 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01043 {
01044 Mm->ip[ipp].stress[idstre] += ipval;
01045 ipp++;
01046 continue;
01047 }
01048 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
01049 {
01050 Mm->ip[ipp].eqother[idoth] += ipval;
01051 ipp++;
01052 continue;
01053 }
01054 if ((ictn[0] & inicond) && (j < nv))
01055 {
01056 if (Mm->ic[ipp] == NULL)
01057 {
01058 Mm->ic[ipp] = new double[nv-j];
01059 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
01060 }
01061 Mm->ic[ipp][idic] += ipval;
01062 ipp++;
01063 continue;
01064 }
01065 ipp++;
01066 }
01067 destrv(gp1); destrv (gp2); destrv (w);
01068 }
01069 }
01070 ipp=Mt->elements[eid].ipp[ri][ci];
01071 ncompstr = Mm->ip[ipp].ncompstr;
01072 ncompeqother = Mm->ip[ipp].ncompeqother;
01073 if ((ictn[0] & inistrain) && (j < ncompstr))
01074 {
01075 nstra++;
01076 idstra++;
01077 continue;
01078 }
01079 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01080 {
01081 nstre++;
01082 idstre++;
01083 continue;
01084 }
01085 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
01086 {
01087 idoth++;
01088 continue;
01089 }
01090 if ((ictn[0] & inicond) && (j < nv))
01091 {
01092 idic++;
01093 continue;
01094 }
01095 }
01096 }