00001 #include "shellq.h"
00002 #include "global.h"
00003 #include "globmat.h"
00004 #include "genfile.h"
00005 #include "node.h"
00006 #include "element.h"
00007 #include "plelemrotlq.h"
00008 #include "q4plate.h"
00009 #include "intpoints.h"
00010 #include <math.h>
00011
00012 shellq::shellq (void)
00013 {
00014 long i,j;
00015
00016 if (Perlq==NULL) Perlq = new planeelemrotlq;
00017 if (Q4pl==NULL) Q4pl = new q4plate;
00018
00019
00020 nne=4;
00021
00022 ndofe=24;
00023
00024 tncomp=6;
00025
00026 napfun=6;
00027
00028 ned=4;
00029
00030 nned=2;
00031
00032 ssst=shell;
00033
00034 nb = Perlq->nb + Q4pl ->nb;
00035
00036 ncomps = new long [Perlq->nb];
00037 ncomps[0]=2;
00038 ncomps[1]=1;
00039
00040 ncompd = new long [Q4pl->nb];
00041 ncompd[0]=3;
00042
00043 nip = new long* [nb];
00044 intordsm = new long* [nb];
00045 for (i=0;i<nb;i++){
00046 nip[i]=new long [nb];
00047 intordsm[i] = new long [nb];
00048 }
00049 tncomps=Perlq->tncomp;
00050 tncompd=Q4pl->tncomp;
00051
00052 long a00,a01,a10,a11,b00,b01,b10,b11;
00053
00054 a00=Perlq->nip[0][0]; a01=Perlq->nip[0][1];
00055 a10=Perlq->nip[1][0]; a11=Perlq->nip[1][1];
00056 b00=Q4pl->nip[0][0]; b01=Q4pl->nip[0][1];
00057 b10=Q4pl->nip[1][0]; b11=Q4pl->nip[1][1];
00058
00059 nip[0][0]=a00; nip[0][1]=a01; nip[0][2]=0; nip[0][3]=0;
00060 nip[1][0]=a10; nip[1][1]=a11; nip[1][2]=0; nip[1][3]=0;
00061 nip[2][0]=0; nip[2][1]=0; nip[2][2]=b00; nip[2][3]=b01;
00062 nip[3][0]=0; nip[3][1]=0; nip[3][2]=b10; nip[3][3]=b11;
00063
00064 a00=Perlq->intordsm[0][0]; a01=Perlq->intordsm[0][1];
00065 a10=Perlq->intordsm[1][0]; a11=Perlq->intordsm[1][1];
00066 b00=Q4pl->intordsm[0][0]; b01=Q4pl->intordsm[0][1];
00067 b10=Q4pl->intordsm[1][0]; b11=Q4pl->intordsm[1][1];
00068
00069 intordsm[0][0]=a00; intordsm[0][1]=a01; intordsm[0][2]=0; intordsm[0][3]=0;
00070 intordsm[1][0]=a10; intordsm[1][1]=a11; intordsm[1][2]=0; intordsm[1][3]=0;
00071 intordsm[2][0]=0; intordsm[2][1]=0; intordsm[2][2]=b00; intordsm[2][3]=b01;
00072 intordsm[3][0]=0; intordsm[3][1]=0; intordsm[3][2]=b10; intordsm[3][3]=b11;
00073
00074 dofe = new long [2];
00075 dofe[0]=Perlq->ndofe;
00076 dofe[1]=Q4pl->ndofe;
00077
00078 tnip=0;
00079 for (i=0;i<nb;i++){
00080 for (j=0;j<nb;j++){
00081 tnip+=nip[i][j];
00082 }
00083 }
00084 }
00085
00086 shellq::~shellq (void)
00087 {
00088 long i;
00089
00090 for (i=0;i<nb;i++){
00091 delete [] nip[i];
00092 delete [] intordsm[i];
00093 }
00094 delete [] intordsm;
00095 delete [] nip;
00096 delete [] ncomps;
00097 delete [] ncompd;
00098 delete [] dofe;
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 void shellq::tran_mat(vector &x, vector &y, matrix &tran, vector &gx, vector &gy, vector &gz)
00123 {
00124 long i,i1;
00125 double dl;
00126 matrix a(3,3);
00127
00128 for (i=0; i<3; i++) {
00129 i1=i+1; if(i1>2)i1=i1-3;
00130
00131 a[0][i]=gx[i1]-gx[i];
00132
00133 a[1][i]=gy[i1]-gy[i];
00134
00135 a[2][i]=gz[i1]-gz[i];
00136 }
00137
00138 dl=sqrt(a[0][0]*a[0][0]+a[1][0]*a[1][0]+a[2][0]*a[2][0]);
00139
00140 tran[0][0]=a[0][0]/dl;
00141 tran[1][0]=a[1][0]/dl;
00142 tran[2][0]=a[2][0]/dl;
00143
00144 tran[0][2]=a[1][0]*a[2][1]-a[2][0]*a[1][1];
00145 tran[1][2]=a[2][0]*a[0][1]-a[0][0]*a[2][1];
00146 tran[2][2]=a[0][0]*a[1][1]-a[1][0]*a[0][1];
00147
00148 dl=sqrt(tran[0][2]*tran[0][2]+tran[1][2]*tran[1][2]+tran[2][2]*tran[2][2]);
00149
00150 tran[0][2]=tran[0][2]/dl;
00151 tran[1][2]=tran[1][2]/dl;
00152 tran[2][2]=tran[2][2]/dl;
00153
00154 tran[0][1]=tran[1][2]*tran[2][0]-tran[2][2]*tran[1][0];
00155 tran[1][1]=tran[2][2]*tran[0][0]-tran[0][2]*tran[2][0];
00156 tran[2][1]=tran[0][2]*tran[1][0]-tran[1][2]*tran[0][0];
00157
00158 dl=sqrt(tran[0][1]*tran[0][1]+tran[1][1]*tran[1][1]+tran[2][1]*tran[2][1]);
00159
00160 tran[0][1]=tran[0][1]/dl;
00161 tran[1][1]=tran[1][1]/dl;
00162 tran[2][1]=tran[2][1]/dl;
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 x[0]=0.0;
00174 y[0]=0.0;
00175 x[1]=sqrt(a[0][0]*a[0][0]+a[1][0]*a[1][0]+a[2][0]*a[2][0]);
00176 y[1]=0.0;
00177 x[2]=(gx(2)-gx(0))*tran(0,0)+(gy(2)-gy(0))*tran(1,0)+(gz(2)-gz(0))*tran(2,0);
00178 y[2]=(gx(2)-gx(0))*tran(0,1)+(gy(2)-gy(0))*tran(1,1)+(gz(2)-gz(0))*tran(2,1);
00179 x[3]=(gx(3)-gx(0))*tran(0,0)+(gy(3)-gy(0))*tran(1,0)+(gz(3)-gz(0))*tran(2,0);
00180 y[3]=(gx(3)-gx(0))*tran(0,1)+(gy(3)-gy(0))*tran(1,1)+(gz(3)-gz(0))*tran(2,1);
00181 }
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194 void shellq::transf_matrix (ivector &nodes,matrix &tmat)
00195 {
00196 long i,i6,n,m;
00197
00198 fillm (0.0,tmat);
00199
00200 n=nodes.n;
00201 m=tmat.m;
00202 for (i=0;i<m;i++){
00203 tmat[i][i]=1.0;
00204 }
00205
00206 for (i=0;i<n;i++){
00207 if (Mt->nodes[nodes[i]].transf>0){
00208 i6=i*6;
00209 tmat[i6][i6] = Mt->nodes[nodes[i]].e1[0]; tmat[i6][i6+1] = Mt->nodes[nodes[i]].e2[0]; tmat[i6][i6+2] = Mt->nodes[nodes[i]].e3[0];
00210 tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1]; tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1]; tmat[i6+1][i6+2] = Mt->nodes[nodes[i]].e3[1];
00211 tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2]; tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2]; tmat[i6+2][i6+2] = Mt->nodes[nodes[i]].e3[2];
00212 i6=i*6+3;
00213 tmat[i6][i6] = Mt->nodes[nodes[i]].e1[0]; tmat[i6][i6+1] = Mt->nodes[nodes[i]].e2[0]; tmat[i6][i6+2] = Mt->nodes[nodes[i]].e3[0];
00214 tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1]; tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1]; tmat[i6+1][i6+2] = Mt->nodes[nodes[i]].e3[1];
00215 tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2]; tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2]; tmat[i6+2][i6+2] = Mt->nodes[nodes[i]].e3[2];
00216 }
00217 }
00218 }
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 void shellq::res_stiffness_matrix (long eid,matrix &sm)
00230 {
00231 long *rcn,*ccn,transf;
00232 vector gx(nne),gy(nne),gz(nne),x(nne),y(nne);
00233 ivector nodes(nne);
00234 matrix lsm,tran(3,3);
00235
00236 Mt->give_elemnodes (eid,nodes);
00237 Mt->give_node_coord3d (gx,gy,gz,eid);
00238
00239 tran_mat(x,y, tran, gx, gy, gz);
00240
00241 fillm (0.0,sm);
00242
00243 rcn = new long [dofe[0]];
00244 ccn = new long [dofe[0]];
00245 allocm (dofe[0],dofe[0],lsm);
00246
00247
00248 Perlq->stiffness_matrix (eid,0,0,lsm,x,y);
00249
00250
00251
00252
00253
00254
00255 rcn[0]=1; rcn[1]=2; rcn[2]=6; rcn[3]=7; rcn[4]=8; rcn[5]=12; rcn[6]=13; rcn[7]=14; rcn[8]=18; rcn[9]=19; rcn[10]=20; rcn[11]=24;
00256 ccn[0]=1; ccn[1]=2; ccn[2]=6; ccn[3]=7; ccn[4]=8; ccn[5]=12; ccn[6]=13; ccn[7]=14; ccn[8]=18; ccn[9]=19; ccn[10]=20; ccn[11]=24;
00257
00258 mat_localize (sm,lsm,rcn,ccn);
00259
00260 delete [] ccn; delete [] rcn; destrm (lsm);
00261
00262 rcn = new long [dofe[1]];
00263 ccn = new long [dofe[1]];
00264 allocm (dofe[1],dofe[1],lsm);
00265
00266
00267 Q4pl->stiffness_matrix (eid,2,2,lsm,x,y);
00268
00269
00270
00271
00272
00273 rcn[0]=3; rcn[1]=4; rcn[2]=5; rcn[3]=9; rcn[4]=10; rcn[5]=11; rcn[6]=15; rcn[7]=16; rcn[8]=17; rcn[9]=21; rcn[10]=22; rcn[11]=23;
00274 ccn[0]=3; ccn[1]=4; ccn[2]=5; ccn[3]=9; ccn[4]=10; ccn[5]=11; ccn[6]=15; ccn[7]=16; ccn[8]=17; ccn[9]=21; ccn[10]=22; ccn[11]=23;
00275
00276 mat_localize (sm,lsm,rcn,ccn);
00277
00278 lgmatrixtransfblock(sm, tran);
00279
00280 transf = Mt->locsystems (nodes);
00281 if (transf>0){
00282 matrix tmat (ndofe,ndofe);
00283 transf_matrix (nodes,tmat);
00284 glmatrixtransf (sm,tmat);
00285 }
00286 delete [] ccn; delete [] rcn; destrm (lsm);
00287
00288 }
00289
00290
00291 void shellq::res_ip_strains (long lcid,long eid)
00292 {
00293 long *rcn,i;
00294 vector aux;
00295 vector gx(nne),gy(nne),gz(nne),x(nne),y(nne),r(ndofe),rs(dofe[0]),rd(dofe[1]);
00296 ivector nodes(nne);
00297 matrix lsm,tmat,tran(3,3);
00298
00299 Mt->give_node_coord3d (gx,gy,gz,eid);
00300 Mt->give_elemnodes (eid,nodes);
00301
00302 rcn = new long [dofe[0]];
00303 allocm (dofe[0],dofe[0],lsm);
00304 eldispl (lcid,eid,r.a);
00305
00306
00307 long transf = Mt->locsystems (nodes);
00308 if (transf>0){
00309 allocv (ndofe,aux);
00310 allocm (ndofe,ndofe,tmat);
00311 transf_matrix (nodes,tmat);
00312
00313 lgvectortransf (aux,r,tmat);
00314 copyv (aux,r);
00315 destrv (aux);
00316 destrm (tmat);
00317 }
00318
00319 tran_mat(x,y, tran, gx, gy, gz);
00320 glvectortransfblock (r, tran);
00321
00322 rcn[0]=1; rcn[1]=2; rcn[2]=6; rcn[3]=7; rcn[4]=8; rcn[5]=12;
00323 rcn[6]=13; rcn[7]=14; rcn[8]=18; rcn[9]=19; rcn[10]=20; rcn[11]=24;
00324 for (i=0;i<Perlq->ndofe;i++){
00325 if (rcn[i]==0) rs[i]=0.0;
00326 else rs[i]=r[rcn[i]-1];
00327 }
00328
00329 Perlq->ip_strains (lcid,eid,0,0,x,y,rs);
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 rcn[0]=3; rcn[1]=4; rcn[2]=5; rcn[3]=9; rcn[4]=10; rcn[5]=11;
00349 rcn[6]=15; rcn[7]=16; rcn[8]=17; rcn[9]=21; rcn[10]=22; rcn[11]=23;
00350 for (i=0;i<Q4pl->ndofe;i++){
00351 if (rcn[i]==0) rd[i]=0.0;
00352 else rd[i]=r[rcn[i]-1];
00353 }
00354
00355 Q4pl->ip_strains (lcid,eid,2,2,x,y,rd);
00356
00357 }
00358
00359
00360 void shellq::nod_strains_ip (long lcid,long eid)
00361 {
00362
00363
00364 Perlq->nod_strains_ip (lcid,eid,0,0);
00365
00366
00367 Q4pl->nod_strains_ip (lcid,eid,2,2);
00368 }
00369
00370 void shellq::strains (long lcid,long eid)
00371 {
00372 }
00373
00374 void shellq::res_ip_stresses (long lcid,long eid)
00375 {
00376 Perlq->compute_nlstress (lcid,eid,0,0);
00377 Q4pl->compute_nlstress (lcid,eid,2,2);
00378 }
00379
00380 void shellq::nod_stresses_ip (long lcid,long eid)
00381 {
00382
00383
00384 Perlq->nod_stresses_ip (lcid,eid,0,0);
00385
00386
00387 Q4pl->nod_stresses_ip (lcid,eid,2,2);
00388 }
00389
00390
00391
00392
00393 void shellq::stresses (long lcid,long eid)
00394 {
00395
00396 }
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422 void shellq::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
00423 {
00424 long i, j, k, l, ipp;
00425 long ii, jj, nv = nodval.n;
00426 double xi, eta;
00427
00428 vector w, gp, anv(nne);
00429 long nstra, nstre, ncompstr, ncompeqother;
00430 long idstra, idstre, idoth, idic;
00431 inictype ict;
00432
00433 nstra = idstra = nstre = idstre = idoth = idic = 0;
00434
00435 ict = ictn[0];
00436 for (i=0; i<nne; i++)
00437 {
00438 if (ictn[i] != ict)
00439 {
00440 print_err("Incompatible types of initial conditions on element %ld\n"
00441 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
00442 abort();
00443 }
00444 }
00445 for (j = 0; j < nv; j++)
00446 {
00447 for(i = 0; i < nne; i++)
00448 anv[i] = nodval[i][j];
00449 for (ii = 0; ii < nb; ii++)
00450 {
00451 for (jj = 0; jj < nb; jj++)
00452 {
00453 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00454 if (intordsm[ii][jj] == 0)
00455 continue;
00456 allocv (intordsm[ii][jj],gp);
00457 allocv (intordsm[ii][jj],w);
00458 gauss_points (gp.a,w.a,intordsm[ii][jj]);
00459 for (k = 0; k < intordsm[ii][jj]; k++)
00460 {
00461 xi=gp[k];
00462 for (l = 0; l < intordsm[ii][jj]; l++)
00463 {
00464 eta=gp[l];
00465
00466
00467 ncompstr = Mm->ip[ipp].ncompstr;
00468 ncompeqother = Mm->ip[ipp].ncompeqother;
00469 if ((ictn[0] & inistrain) && (j < ncompstr))
00470 {
00471
00472 ipp++;
00473 continue;
00474 }
00475 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
00476 {
00477
00478 ipp++;
00479 continue;
00480 }
00481 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
00482 {
00483
00484 ipp++;
00485 continue;
00486 }
00487 if ((ictn[0] & inicond) && (j < nv))
00488 {
00489 if (Mm->ic[ipp] == NULL)
00490 {
00491 Mm->ic[ipp] = new double[nv-j];
00492 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
00493 }
00494
00495 ipp++;
00496 continue;
00497 }
00498 ipp++;
00499 }
00500 }
00501 destrv (gp); destrv (w);
00502 }
00503 }
00504 ipp=Mt->elements[eid].ipp[ri][ci];
00505 ncompstr = Mm->ip[ipp].ncompstr;
00506 ncompeqother = Mm->ip[ipp].ncompeqother;
00507 if ((ictn[0] & inistrain) && (j < ncompstr))
00508 {
00509 nstra++;
00510 idstra++;
00511 continue;
00512 }
00513 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
00514 {
00515 nstre++;
00516 idstre++;
00517 continue;
00518 }
00519 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
00520 {
00521 idoth++;
00522 continue;
00523 }
00524 if ((ictn[0] & inicond) && (j < nv))
00525 {
00526 idic++;
00527 continue;
00528 }
00529 }
00530 }