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