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