00001 #include "plelemqt.h"
00002 #include "global.h"
00003 #include "globmat.h"
00004 #include "genfile.h"
00005 #include "adaptivity.h"
00006 #include "node.h"
00007 #include "element.h"
00008 #include "intpoints.h"
00009 #include "plelemlt.h"
00010 #include "plelemsubqt.h"
00011 #include "plelemlq.h"
00012 #include "plelemqq.h"
00013 #include "loadcase.h"
00014 #include "gadaptivity.h"
00015 #include <stdlib.h>
00016 #include <math.h>
00017
00018
00019 planeelemqt::planeelemqt (void)
00020 {
00021 long i,j;
00022
00023
00024 nne=6;
00025
00026 ndofe=12;
00027
00028 tncomp=3;
00029
00030 gncomp=4;
00031
00032 napfun=2;
00033
00034 intordmm=6;
00035
00036 ned=3;
00037
00038 nned=3;
00039
00040 intordb=3;
00041
00042
00043 nb=2;
00044
00045
00046 ncomp = new long [nb];
00047 ncomp[0]=2;
00048 ncomp[1]=1;
00049
00050
00051 cncomp = new long [nb];
00052 cncomp[0]=0;
00053 cncomp[1]=2;
00054
00055
00056
00057 nip = new long* [nb];
00058 intordsm = new long* [nb];
00059 for (i=0;i<nb;i++){
00060 nip[i] = new long [nb];
00061 intordsm[i] = new long [nb];
00062 }
00063
00064 nip[0][0]=3; nip[0][1]=0;
00065 nip[1][0]=0; nip[1][1]=3;
00066
00067
00068 tnip=0;
00069 for (i=0;i<nb;i++){
00070 for (j=0;j<nb;j++){
00071 tnip+=nip[i][j];
00072 }
00073 }
00074
00075 intordsm[0][0]=3; intordsm[0][1]=0;
00076 intordsm[1][0]=0; intordsm[1][1]=3;
00077
00078 }
00079
00080 planeelemqt::~planeelemqt (void)
00081 {
00082 long i;
00083
00084 for (i=0;i<nb;i++){
00085 delete [] nip[i];
00086 delete [] intordsm[i];
00087 }
00088 delete [] nip;
00089 delete [] intordsm;
00090
00091 delete [] ncomp;
00092 delete [] cncomp;
00093 }
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 double planeelemqt::approx (double xi,double eta,vector &nodval)
00106 {
00107 double f;
00108 vector bf(nne);
00109
00110 bf_quad_3_2d (bf.a,xi,eta);
00111 scprd (bf,nodval,f);
00112 return f;
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123 void planeelemqt::bf_matrix (matrix &n,double xi,double eta)
00124 {
00125 long i,i1,i2;
00126 vector bf(nne);
00127
00128 bf_quad_3_2d (bf.a,xi,eta);
00129 fillm (0.0,n);
00130
00131 i1=0; i2=1;
00132 for (i=0;i<nne;i++){
00133 n[0][i1]=bf[i]; i1+=2;
00134 n[1][i2]=bf[i]; i2+=2;
00135 }
00136 }
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 void planeelemqt::geom_matrix (matrix &gm,vector &x,vector &y,double xi,double eta,double &jac)
00149 {
00150 long i,i1,i2;
00151 vector dx(nne),dy(nne);
00152
00153 dx_bf_quad_3_2d (dx.a,xi,eta);
00154 dy_bf_quad_3_2d (dy.a,xi,eta);
00155
00156 derivatives_2d (dx,dy,jac,x,y,xi,eta);
00157
00158 fillm (0.0,gm);
00159
00160 i1=0; i2=1;
00161 for (i=0;i<nne;i++){
00162 gm[0][i1]=dx[i];
00163 gm[1][i2]=dy[i];
00164 gm[2][i1]=dy[i]; i1+=2;
00165 gm[2][i2]=dx[i]; i2+=2;
00166 }
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 void planeelemqt::geom_matrix_block (matrix &gm,long ri,vector &x,vector &y,double xi,double eta,double &jac)
00180 {
00181 long i,i1,i2;
00182 vector dx(nne),dy(nne);
00183
00184 dx_bf_quad_3_2d (dx.a,xi,eta);
00185 dy_bf_quad_3_2d (dy.a,xi,eta);
00186
00187 derivatives_2d (dx,dy,jac,x,y,xi,eta);
00188
00189 fillm (0.0,gm);
00190
00191 if (ri==0){
00192 i1=0; i2=1;
00193 for (i=0;i<nne;i++){
00194 gm[0][i1]=dx[i]; i1+=2;
00195 gm[1][i2]=dy[i]; i2+=2;
00196 }
00197 }
00198
00199 if (ri==1){
00200 i1=0; i2=1;
00201 for (i=0;i<nne;i++){
00202 gm[0][i1]=dy[i]; i1+=2;
00203 gm[0][i2]=dx[i]; i2+=2;
00204 }
00205 }
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 void planeelemqt::dmatblock (long ri,long ci,matrix &d, matrix &dd)
00217 {
00218 fillm (0.0,dd);
00219
00220 if (ri==0 && ci==0){
00221 dd[0][0]=d[0][0]; dd[0][1]=d[0][1];
00222 dd[1][0]=d[1][0]; dd[1][1]=d[1][1];
00223 }
00224 if (ri==0 && ci==1){
00225 dd[0][0]=d[0][2];
00226 dd[1][0]=d[1][2];
00227 }
00228 if (ri==1 && ci==0){
00229 dd[0][0]=d[2][0]; dd[0][1]=d[2][1];
00230 }
00231 if (ri==1 && ci==1){
00232 dd[0][0]=d[2][2];
00233 }
00234 }
00235
00236
00237
00238
00239
00240
00241 void planeelemqt::transf_matrix (ivector &nodes,matrix &tmat)
00242 {
00243 long i,n,m;
00244
00245 fillm (0.0,tmat);
00246
00247 n=nodes.n;
00248 m=tmat.m;
00249 for (i=0;i<m;i++){
00250 tmat[i][i]=1.0;
00251 }
00252
00253 for (i=0;i<n;i++){
00254 if (Mt->nodes[nodes[i]].transf>0){
00255 tmat[i*2][i*2] = Mt->nodes[nodes[i]].e1[0]; tmat[i*2][i*2+1] = Mt->nodes[nodes[i]].e2[0];
00256 tmat[i*2+1][i*2] = Mt->nodes[nodes[i]].e1[1]; tmat[i*2+1][i*2+1] = Mt->nodes[nodes[i]].e2[1];
00257 }
00258 }
00259 }
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270 void planeelemqt::stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
00271 {
00272 long i,ii,jj,ipp;
00273 double jac,thick;
00274 ivector nodes(nne);
00275 vector t(nne),gp1,gp2,w;
00276 matrix gmr,gmc,dd,d(tncomp,tncomp);
00277
00278 Mt->give_elemnodes (eid,nodes);
00279 Mc->give_thickness (eid,nodes,t);
00280
00281 fillm (0.0,sm);
00282
00283 for (ii=0;ii<nb;ii++){
00284 allocm (ncomp[ii],ndofe,gmr);
00285 for (jj=0;jj<nb;jj++){
00286 if (intordsm[ii][jj]==0) continue;
00287
00288 allocm (ncomp[jj],ndofe,gmc);
00289 allocm (ncomp[ii],ncomp[jj],dd);
00290 allocv (intordsm[ii][jj],gp1);
00291 allocv (intordsm[ii][jj],gp2);
00292 allocv (intordsm[ii][jj],w);
00293
00294 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
00295 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00296
00297 for (i=0;i<intordsm[ii][jj];i++){
00298
00299 geom_matrix_block (gmr,ii,x,y,gp1[i],gp2[i],jac);
00300 geom_matrix_block (gmc,jj,x,y,gp1[i],gp2[i],jac);
00301
00302
00303 Mm->matstiff (d,ipp);
00304 dmatblock (ii,jj,d,dd);
00305
00306
00307 thick = approx (gp1[i],gp2[i],t);
00308
00309 jac*=w[i]*thick;
00310
00311
00312
00313
00314
00315 bdbjac (sm,gmr,dd,gmc,jac);
00316
00317 ipp++;
00318 }
00319 destrm (dd); destrm (gmc); destrv (gp1); destrv (gp2); destrv (w);
00320 }
00321 destrm (gmr);
00322 }
00323
00324 }
00325
00326 void planeelemqt::res_stiffness_matrix (long eid,matrix &sm)
00327 {
00328 long transf;
00329 vector x(nne),y(nne);
00330 ivector nodes(nne);
00331 Mt->give_elemnodes (eid,nodes);
00332 Mt->give_node_coord2d (x,y,eid);
00333
00334 stiffness_matrix (eid,0,0,sm,x,y);
00335
00336
00337 transf = Mt->locsystems (nodes);
00338 if (transf>0){
00339 matrix tmat (ndofe,ndofe);
00340 transf_matrix (nodes,tmat);
00341 glmatrixtransf (sm,tmat);
00342 }
00343 }
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 void planeelemqt::mass_matrix (long eid,matrix &mm)
00355 {
00356 long i;
00357 double jac,thick,rho;
00358 ivector nodes(nne);
00359 vector x(nne),y(nne),w(intordmm),gp1(intordmm),gp2(intordmm),t(nne),dens(nne);
00360 matrix n(napfun,ndofe);
00361
00362 Mt->give_elemnodes (eid,nodes);
00363 Mt->give_node_coord2d (x,y,eid);
00364 Mc->give_thickness (eid,nodes,t);
00365 Mc->give_density (eid,nodes,dens);
00366
00367 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00368
00369 fillm (0.0,mm);
00370
00371 for (i=0;i<intordmm;i++){
00372
00373 bf_matrix (n,gp1[i],gp2[i]);
00374
00375
00376 thick = approx (gp1[i],gp2[i],t);
00377
00378
00379 rho = approx (gp1[i],gp2[i],dens);
00380
00381
00382 jac_2d (jac,x,y,gp1[i],gp2[i]);
00383
00384 jac*=w[i]*thick*rho;
00385
00386
00387 nnj (mm.a,n.a,jac,n.m,n.n);
00388 }
00389
00390 }
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 void planeelemqt::load_matrix (long eid,matrix &lm)
00404 {
00405 long i;
00406 double jac,thick;
00407 ivector nodes(nne);
00408 vector x(nne),y(nne),w(intordmm),gp1(intordmm),gp2(intordmm),t(nne);
00409 matrix n(napfun,ndofe);
00410
00411 Mt->give_elemnodes (eid,nodes);
00412 Mt->give_node_coord2d (x,y,eid);
00413 Mc->give_thickness (eid,nodes,t);
00414
00415 gauss_points_tr (gp1.a,gp2.a,w.a,intordmm);
00416
00417 fillm (0.0,lm);
00418
00419 for (i=0;i<intordmm;i++){
00420 bf_matrix (n,gp1[i],gp2[i]);
00421
00422 thick = approx (gp1[i],gp2[i],t);
00423
00424 jac_2d (jac,x,y,gp1[i],gp2[i]);
00425
00426 jac*=w[i]*thick;
00427
00428 nnj (lm.a,n.a,jac,n.m,n.n);
00429 }
00430
00431 }
00432
00433
00434
00435
00436
00437
00438 void planeelemqt::res_mainip_strains (long lcid,long eid)
00439 {
00440 vector aux,x(nne),y(nne),r(ndofe);
00441 ivector nodes(nne);
00442 matrix tmat;
00443
00444 Mt->give_node_coord2d (x,y,eid);
00445 Mt->give_elemnodes (eid,nodes);
00446 eldispl (lcid,eid,r.a);
00447
00448
00449 long transf = Mt->locsystems (nodes);
00450 if (transf>0){
00451 allocv (ndofe,aux);
00452 allocm (ndofe,ndofe,tmat);
00453 transf_matrix (nodes,tmat);
00454
00455 lgvectortransf (aux,r,tmat);
00456 copyv (aux,r);
00457 destrv (aux);
00458 destrm (tmat);
00459 }
00460
00461 mainip_strains (lcid,eid,0,0,x,y,r);
00462 }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474 void planeelemqt::mainip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
00475 {
00476 long i,ii,ipp;
00477 double jac;
00478 vector gp1,gp2,w,eps;
00479 matrix gm;
00480
00481 for (ii=0;ii<nb;ii++){
00482 allocv (intordsm[ii][ii],gp1);
00483 allocv (intordsm[ii][ii],gp2);
00484 allocv (intordsm[ii][ii],w);
00485 allocv (ncomp[ii],eps);
00486 allocm (ncomp[ii],ndofe,gm);
00487
00488 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00489
00490 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00491 for (i=0;i<intordsm[ii][ii];i++){
00492 geom_matrix_block (gm,ii,x,y,gp1[i],gp2[i],jac);
00493 mxv (gm,r,eps);
00494
00495 Mm->storestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00496 ipp++;
00497 }
00498 destrm (gm); destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
00499 }
00500 }
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510 void planeelemqt::nod_strains (long lcid,long eid,long ri,long ci)
00511 {
00512 long i,ii,ipp;
00513 double *lsm,*lhs,*rhs;
00514 vector nxi(nne),neta(nne),gp1,gp2,w,eps,aux,natcoord(2);
00515 ivector nodes(nne);
00516
00517 lsm = new double [9];
00518
00519 nodecoord (nxi,neta);
00520 Mt->give_elemnodes (eid,nodes);
00521
00522 for (ii=0;ii<nb;ii++){
00523 allocv (intordsm[ii][ii],gp1);
00524 allocv (intordsm[ii][ii],gp2);
00525 allocv (intordsm[ii][ii],w);
00526 allocv (ncomp[ii],eps);
00527 lhs = new double [ncomp[ii]*3];
00528 rhs = new double [ncomp[ii]*3];
00529 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00530
00531 nullv (lsm,9);
00532 nullv (rhs,ncomp[ii]*3);
00533
00534 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00535 for (i=0;i<intordsm[ii][ii];i++){
00536 Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00537
00538 natcoord[0]=gp1[i]; natcoord[1]=gp2[i];
00539 matassem_lsm (lsm,natcoord);
00540 rhsassem_lsm (rhs,natcoord,eps);
00541
00542 ipp++;
00543 }
00544
00545 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
00546 Mt->strain_nodal_values (nodes,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii],lcid);
00547
00548 delete [] lhs; delete [] rhs;
00549 destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
00550 }
00551 delete [] lsm;
00552 }
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562 void planeelemqt::elem_strains (double **stra,long lcid,long eid,long ri,long ci)
00563 {
00564 long i,ii,ipp;
00565 double *lsm,*lhs,*rhs;
00566 vector nxi(nne),neta(nne),gp1,gp2,w,eps,aux,natcoord(2);
00567
00568 lsm = new double [9];
00569
00570 nodecoord (nxi,neta);
00571
00572 for (ii=0;ii<nb;ii++){
00573 allocv (intordsm[ii][ii],gp1);
00574 allocv (intordsm[ii][ii],gp2);
00575 allocv (intordsm[ii][ii],w);
00576 allocv (ncomp[ii],eps);
00577 lhs = new double [ncomp[ii]*3];
00578 rhs = new double [ncomp[ii]*3];
00579 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00580
00581 nullv (lsm,9);
00582 nullv (rhs,ncomp[ii]*3);
00583
00584 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00585 for (i=0;i<intordsm[ii][ii];i++){
00586 Mm->givestrain (lcid,ipp,cncomp[ii],ncomp[ii],eps);
00587
00588 natcoord[0]=gp1[i]; natcoord[1]=gp2[i];
00589 matassem_lsm (lsm,natcoord);
00590 rhsassem_lsm (rhs,natcoord,eps);
00591
00592 ipp++;
00593 }
00594
00595 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
00596 nodal_values (stra,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii]);
00597
00598 delete [] lhs; delete [] rhs;
00599 destrv (eps); destrv (w); destrv (gp1); destrv (gp2);
00600 }
00601
00602 delete [] lsm;
00603 }
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 void planeelemqt::appstrain (long lcid,long eid,double xi,double eta,long fi,long ncomp,vector &eps)
00617 {
00618 long i,j,k;
00619 ivector nodes;
00620 vector nodval;
00621
00622 if (ncomp != eps.n){
00623 fprintf (stderr,"\n\n wrong interval of indices in function strain (%s, line %d).\n",__FILE__,__LINE__);
00624 abort ();
00625 }
00626
00627 allocv (nne,nodes);
00628 allocv (nne,nodval);
00629 Mt->give_elemnodes (eid,nodes);
00630 k=0;
00631 for (i=fi;i<fi+ncomp;i++){
00632 for (j=0;j<nne;j++){
00633 nodval[j]=Mt->nodes[nodes[j]].strain[lcid*tncomp+i];
00634 }
00635 eps[k]=approx (xi,eta,nodval);
00636 k++;
00637 }
00638
00639 destrv (nodes); destrv (nodval);
00640 }
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651 void planeelemqt::allip_strains (double **stra,long lcid,long eid,long ri,long ci)
00652 {
00653 long i,ii,jj,ipp;
00654 vector eps,gp1,gp2,w;
00655
00656 allocv (tncomp,eps);
00657
00658 for (ii=0;ii<nb;ii++){
00659 for (jj=0;jj<nb;jj++){
00660 if (intordsm[ii][jj]==0) continue;
00661
00662 allocv (intordsm[ii][jj],gp1);
00663 allocv (intordsm[ii][jj],gp2);
00664 allocv (intordsm[ii][jj],w);
00665 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
00666 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
00667
00668 for (i=0;i<intordsm[ii][jj];i++){
00669
00670 if (Mp->strainaver==0)
00671 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
00672 if (Mp->strainaver==1)
00673 appstrain (lcid,eid,gp1[i],gp2[i],0,tncomp,eps);
00674
00675 Mm->storestrain (lcid,ipp,eps);
00676 ipp++;
00677 }
00678 destrv (w); destrv (gp2); destrv (gp1);
00679 }
00680 }
00681 destrv (eps);
00682 }
00683
00684 void planeelemqt::strains (long lcid,long eid,long ri,long ci)
00685 {
00686 long i,naep,ncp,sid;
00687 double **stra;
00688 vector coord,eps;
00689
00690 if (Mp->strainaver==0){
00691 stra = new double* [nne];
00692 for (i=0;i<nne;i++){
00693 stra[i] = new double [tncomp];
00694 }
00695 elem_strains (stra,lcid,eid,ri,ci);
00696 }
00697
00698 switch (Mm->stra.tape[eid]){
00699 case nowhere:{
00700 break;
00701 }
00702 case intpts:{
00703 allip_strains (stra,lcid,eid,ri,ci);
00704 break;
00705 }
00706 case enodes:{
00707 break;
00708 }
00709 case userdefined:{
00710
00711 naep = Mm->stra.give_naep (eid);
00712 ncp = Mm->stra.give_ncomp (eid);
00713 sid = Mm->stra.give_sid (eid);
00714 allocv (ncp,eps);
00715 allocv (2,coord);
00716 for (i=0;i<naep;i++){
00717 Mm->stra.give_aepcoord (sid,i,coord);
00718
00719 if (Mp->strainaver==0)
00720 appval (coord[0],coord[1],0,ncp,eps,stra);
00721 if (Mp->strainaver==1)
00722 appstrain (lcid,eid,coord[0],coord[1],0,ncp,eps);
00723
00724 Mm->stra.storevalues(lcid,eid,i,eps);
00725 }
00726 destrv (eps);
00727 destrv (coord);
00728 break;
00729 }
00730 default:{
00731 fprintf (stderr,"\n\n unknown strain point is required in function planeelemlt::strains (%s, line %d).\n",__FILE__,__LINE__);
00732 }
00733 }
00734
00735 if (Mp->strainaver==0){
00736 for (i=0;i<nne;i++){
00737 delete [] stra[i];
00738 }
00739 delete [] stra;
00740 }
00741 }
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754 void planeelemqt::nodecoord (vector &xi,vector &eta)
00755 {
00756 xi[0] = 1.0; eta[0] = 0.0;
00757 xi[1] = 0.0; eta[1] = 1.0;
00758 xi[2] = 0.0; eta[2] = 0.0;
00759 xi[3] = 0.5; eta[3] = 0.5;
00760 xi[4] = 0.0; eta[4] = 0.5;
00761 xi[5] = 0.5; eta[5] = 0.0;
00762 }
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773 void planeelemqt::appval (double xi,double eta,long fi,long nc,vector &eps,double **val)
00774 {
00775 long i,j,k;
00776 vector nodval;
00777
00778 k=0;
00779 allocv (nne,nodval);
00780 for (i=fi;i<fi+nc;i++){
00781 for (j=0;j<nne;j++){
00782 nodval[j]=val[j][i];
00783 }
00784 eps[k]=approx (xi,eta,nodval);
00785 k++;
00786 }
00787
00788 destrv (nodval);
00789 }
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802 void planeelemqt::mainip_stresses (long lcid,long eid,long ri,long ci)
00803 {
00804 long i,ii,jj,ipp;
00805 vector gp1,gp2,w,eps,epst,epstt,sig,auxsig;
00806 matrix gm,d(tncomp,tncomp),dd;
00807
00808 allocm (tncomp,tncomp,d);
00809
00810 for (ii=0;ii<nb;ii++){
00811 allocv (ncomp[ii],sig);
00812 allocv (ncomp[ii],auxsig);
00813 allocv (intordsm[ii][ii],gp1);
00814 allocv (intordsm[ii][ii],gp2);
00815 allocv (intordsm[ii][ii],w);
00816
00817 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00818 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00819
00820 for (i=0;i<intordsm[ii][ii];i++){
00821 Mm->matstiff (d,ipp);
00822
00823 fillv (0.0,sig);
00824 for (jj=0;jj<nb;jj++){
00825 allocv (ncomp[jj],eps);
00826 allocm (ncomp[ii],ncomp[jj],dd);
00827
00828 if (Mp->strainaver==0)
00829 Mm->givestrain (lcid,ipp,cncomp[jj],ncomp[jj],eps);
00830 if (Mp->strainaver==1)
00831 appstrain (lcid,eid,gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps);
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 dmatblock (ii,jj,d,dd);
00845 mxv (dd,eps,auxsig);
00846 addv (auxsig,sig,sig);
00847 destrm (dd); destrv (eps);
00848 }
00849
00850 Mm->storestress (lcid,ipp,sig);
00851
00852 ipp++;
00853 }
00854 destrv (w); destrv (gp2); destrv (gp1); destrv (auxsig); destrv (sig);
00855 }
00856
00857 }
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868 void planeelemqt::nod_stresses (long lcid,long eid,long ri,long ci)
00869 {
00870 long i,ii,jj,ipp;
00871 double *lsm,*lhs,*rhs;
00872 vector nxi(nne),neta(nne),gp1,gp2,w,eps,epst,epstt,sig,auxsig,natcoord(2);
00873 ivector nodes(nne);
00874 matrix d(tncomp,tncomp),dd;
00875
00876 lsm = new double [9];
00877
00878 Mt->give_elemnodes (eid,nodes);
00879
00880 nodecoord (nxi,neta);
00881
00882 for (ii=0;ii<nb;ii++){
00883 allocv (intordsm[ii][ii],gp1);
00884 allocv (intordsm[ii][ii],gp2);
00885 allocv (intordsm[ii][ii],w);
00886 allocv (ncomp[ii],sig);
00887 allocv (ncomp[ii],auxsig);
00888 lhs = new double [ncomp[ii]*3];
00889 rhs = new double [ncomp[ii]*3];
00890 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00891
00892 nullv (lsm,9);
00893 nullv (rhs,ncomp[ii]*3);
00894
00895 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00896
00897 for (i=0;i<intordsm[ii][ii];i++){
00898
00899 Mm->matstiff (d,ipp);
00900
00901 fillv (0.0,sig);
00902 for (jj=0;jj<nb;jj++){
00903 allocv (ncomp[jj],eps);
00904 allocm (ncomp[ii],ncomp[jj],dd);
00905
00906 if (Mp->strainaver==0)
00907 Mm->givestrain (lcid,ipp,cncomp[jj],ncomp[jj],eps);
00908 if (Mp->strainaver==1)
00909 appstrain (lcid,eid,gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps);
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922 dmatblock (ii,jj,d,dd);
00923 mxv (dd,eps,auxsig);
00924 addv (auxsig,sig,sig);
00925 destrm (dd); destrv (eps);
00926 }
00927
00928 natcoord[0]=gp1[i]; natcoord[1]=gp2[i];
00929 matassem_lsm (lsm,natcoord);
00930 rhsassem_lsm (rhs,natcoord,sig);
00931
00932 ipp++;
00933 }
00934
00935 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
00936 Mt->stress_nodal_values (nodes,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii],lcid);
00937
00938
00939 delete [] lhs; delete [] rhs;
00940 destrv (auxsig); destrv (sig); destrv (eps); destrv (w); destrv (gp2); destrv (gp1);
00941 }
00942
00943 delete [] lsm;
00944 }
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955 void planeelemqt::elem_stresses (double **stra,double **stre,long lcid,long eid,long ri,long ci)
00956 {
00957 long i,ii,jj,ipp;
00958 double *lsm,*lhs,*rhs;
00959 vector nxi(nne),neta(nne),gp1,gp2,w,eps,epst,epstt,sig,auxsig,natcoord(2);
00960 matrix d(tncomp,tncomp),dd;
00961
00962 lsm = new double [9];
00963
00964 nodecoord (nxi,neta);
00965
00966 for (ii=0;ii<nb;ii++){
00967 allocv (intordsm[ii][ii],gp1);
00968 allocv (intordsm[ii][ii],gp2);
00969 allocv (intordsm[ii][ii],w);
00970 allocv (ncomp[ii],sig);
00971 allocv (ncomp[ii],auxsig);
00972 lhs = new double [ncomp[ii]*3];
00973 rhs = new double [ncomp[ii]*3];
00974 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
00975
00976 nullv (lsm,9);
00977 nullv (rhs,ncomp[ii]*3);
00978
00979 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
00980
00981 for (i=0;i<intordsm[ii][ii];i++){
00982
00983 Mm->matstiff (d,ipp);
00984
00985 fillv (0.0,sig);
00986 for (jj=0;jj<nb;jj++){
00987 allocv (ncomp[jj],eps);
00988 allocm (ncomp[ii],ncomp[jj],dd);
00989
00990 if (Mp->strainaver==0)
00991 appval (gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps,stra);
00992 if (Mp->strainaver==1)
00993 appstrain (lcid,eid,gp1[i],gp2[i],cncomp[jj],ncomp[jj],eps);
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006 dmatblock (ii,jj,d,dd);
01007 mxv (dd,eps,auxsig);
01008 addv (auxsig,sig,sig);
01009 destrm (dd); destrv (eps);
01010 }
01011
01012 natcoord[0]=gp1[i]; natcoord[1]=gp2[i];
01013 matassem_lsm (lsm,natcoord);
01014 rhsassem_lsm (rhs,natcoord,sig);
01015
01016 ipp++;
01017 }
01018
01019 solve_lsm (lsm,lhs,rhs,Mp->zero,3,ncomp[ii]);
01020 nodal_values (stre,nxi,neta,nxi,lhs,2,cncomp[ii],ncomp[ii]);
01021
01022 delete [] lhs; delete [] rhs;
01023 destrv (auxsig); destrv (sig); destrv (eps); destrv (w); destrv (gp2); destrv (gp1);
01024 }
01025
01026 delete [] lsm;
01027 }
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040 void planeelemqt::appstress (long lcid,long eid,double xi,double eta,long fi,long ncomp,vector &sig)
01041 {
01042 long i,j,k;
01043 ivector nodes;
01044 vector nodval;
01045
01046 if (ncomp != sig.n){
01047 fprintf (stderr,"\n\n wrong interval of indices in function stress (%s, line %d).\n",__FILE__,__LINE__);
01048 abort ();
01049 }
01050
01051 allocv (nne,nodes);
01052 allocv (nne,nodval);
01053
01054 Mt->give_elemnodes (eid,nodes);
01055 k=0;
01056 for (i=fi;i<fi+ncomp;i++){
01057 for (j=0;j<nne;j++){
01058 nodval[j]=Mt->nodes[nodes[j]].stress[lcid*tncomp+i];
01059 }
01060 sig[k]=approx (xi,eta,nodval);
01061 k++;
01062 }
01063
01064 destrv (nodes); destrv (nodval);
01065 }
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076 void planeelemqt::allip_stresses (double **stre,long lcid,long eid,long ri,long ci)
01077 {
01078 long i,ii,jj,ipp;
01079 vector sig(tncomp),gp1,gp2,w;
01080
01081 for (ii=0;ii<nb;ii++){
01082 for (jj=0;jj<nb;jj++){
01083 if (intordsm[ii][jj]==0) continue;
01084
01085 allocv (intordsm[ii][jj],gp1);
01086 allocv (intordsm[ii][jj],gp2);
01087 allocv (intordsm[ii][jj],w);
01088 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
01089 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01090
01091 for (i=0;i<intordsm[ii][jj];i++){
01092
01093 if (Mp->stressaver==0)
01094 appval (gp1[i],gp2[i],0,tncomp,sig,stre);
01095 if (Mp->stressaver==1)
01096 appstress (lcid,eid,gp1[i],gp2[i],0,tncomp,sig);
01097
01098 Mm->storestress (lcid,ipp,sig);
01099 ipp++;
01100 }
01101 destrv (w); destrv (gp2); destrv (gp1);
01102 }
01103 }
01104 }
01105
01106 void planeelemqt::stresses (long lcid,long eid,long ri,long ci)
01107 {
01108 long i,naep,ncp,sid;
01109 double **stra,**stre;
01110 vector coord,sig;
01111
01112 if (Mp->stressaver==0){
01113 stra = new double* [nne];
01114 stre = new double* [nne];
01115 for (i=0;i<nne;i++){
01116 stra[i] = new double [tncomp];
01117 stre[i] = new double [tncomp];
01118 }
01119 elem_strains (stra,lcid,eid,ri,ci);
01120 elem_stresses (stra,stre,lcid,eid,ri,ci);
01121 }
01122
01123 switch (Mm->stre.tape[eid]){
01124 case nowhere:{
01125 break;
01126 }
01127 case intpts:{
01128 allip_stresses (stre,lcid,eid,ri,ci);
01129 break;
01130 }
01131 case enodes:{
01132 break;
01133 }
01134 case userdefined:{
01135
01136 naep = Mm->stre.give_naep (eid);
01137 ncp = Mm->stre.give_ncomp (eid);
01138 sid = Mm->stre.give_sid (eid);
01139 allocv (ncp,sig);
01140 allocv (2,coord);
01141 for (i=0;i<naep;i++){
01142 Mm->stre.give_aepcoord (sid,i,coord);
01143
01144 if (Mp->stressaver==0)
01145 appval (coord[0],coord[1],0,ncp,sig,stre);
01146 if (Mp->stressaver==1)
01147 appstress (lcid,eid,coord[0],coord[1],0,ncp,sig);
01148
01149 Mm->stre.storevalues(lcid,eid,i,sig);
01150 }
01151 destrv (sig);
01152 destrv (coord);
01153 break;
01154 }
01155 default:{
01156 fprintf (stderr,"\n\n unknown stress point is required in function planeelemlq::stresses (%s, line %d).\n",__FILE__,__LINE__);
01157 }
01158 }
01159
01160 if (Mp->stressaver==0){
01161 for (i=0;i<nne;i++){
01162 delete [] stra[i];
01163 delete [] stre[i];
01164 }
01165 delete [] stra;
01166 delete [] stre;
01167 }
01168
01169 }
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190 void planeelemqt::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01191 {
01192 integratedquant iq;
01193 iq=locstress;
01194
01195
01196 compute_nlstress (lcid,eid,ri,ci);
01197
01198 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01199 }
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214 void planeelemqt::nonloc_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01215 {
01216 integratedquant iq;
01217 iq=nonlocstress;
01218
01219
01220 compute_nonloc_nlstress (lcid,eid,ri,ci);
01221
01222 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01223 }
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238 void planeelemqt::incr_internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
01239 {
01240 integratedquant iq;
01241 iq=stressincr;
01242
01243
01244 compute_nlstressincr (lcid,eid,ri,ci);
01245
01246 elem_integration (iq,lcid,eid,ri,ci,ifor,x,y);
01247 }
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261 void planeelemqt::eigstrain_forces (long lcid,long eid,long ri,long ci,vector &nfor,vector &x,vector &y)
01262 {
01263 integratedquant iq;
01264 iq=eigstress;
01265
01266
01267 if (Mp->eigstrcomp)
01268 compute_eigstress (lcid,eid,ri,ci);
01269
01270
01271 elem_integration (iq,lcid,eid,ri,ci,nfor,x,y);
01272 }
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285 void planeelemqt::res_internal_forces (long lcid,long eid,vector &ifor)
01286 {
01287 long transf;
01288 ivector nodes (nne);
01289 vector v(ndofe),x(nne),y(nne);
01290
01291 Mt->give_node_coord2d (x,y,eid);
01292 internal_forces (lcid,eid,0,0,ifor,x,y);
01293
01294
01295 Mt->give_elemnodes (eid,nodes);
01296 transf = Mt->locsystems (nodes);
01297 if (transf>0){
01298 matrix tmat (ndofe,ndofe);
01299 transf_matrix (nodes,tmat);
01300 glvectortransf (ifor,v,tmat);
01301 copyv (v,ifor);
01302 }
01303 }
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316 void planeelemqt::res_nonloc_internal_forces (long lcid,long eid,vector &ifor)
01317 {
01318 long transf;
01319 ivector nodes (nne);
01320 vector v(ndofe),x(nne),y(nne);
01321
01322 Mt->give_node_coord2d (x,y,eid);
01323 nonloc_internal_forces (lcid,eid,0,0,ifor,x,y);
01324
01325
01326 Mt->give_elemnodes (eid,nodes);
01327 transf = Mt->locsystems (nodes);
01328 if (transf>0){
01329 matrix tmat (ndofe,ndofe);
01330 transf_matrix (nodes,tmat);
01331 glvectortransf (ifor,v,tmat);
01332 copyv (v,ifor);
01333 }
01334 }
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347 void planeelemqt::res_incr_internal_forces (long lcid,long eid,vector &ifor)
01348 {
01349 long transf;
01350 ivector nodes (nne);
01351 vector v(ndofe),x(nne),y(nne);
01352
01353 Mt->give_node_coord2d (x,y,eid);
01354 incr_internal_forces (lcid,eid,0,0,ifor,x,y);
01355
01356
01357 Mt->give_elemnodes (eid,nodes);
01358 transf = Mt->locsystems (nodes);
01359 if (transf>0){
01360 matrix tmat (ndofe,ndofe);
01361 transf_matrix (nodes,tmat);
01362 glvectortransf (ifor,v,tmat);
01363 copyv (v,ifor);
01364 }
01365 }
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378 void planeelemqt::res_eigstrain_forces (long lcid,long eid,vector &nfor)
01379 {
01380 long transf;
01381 ivector nodes (nne);
01382 vector v(ndofe),x(nne),y(nne);
01383
01384 Mt->give_node_coord2d (x,y,eid);
01385 eigstrain_forces (lcid,eid,0,0,nfor,x,y);
01386
01387
01388 Mt->give_elemnodes (eid,nodes);
01389 transf = Mt->locsystems (nodes);
01390 if (transf>0){
01391 matrix tmat (ndofe,ndofe);
01392 transf_matrix (nodes,tmat);
01393 glvectortransf (nfor,v,tmat);
01394 copyv (v,nfor);
01395 }
01396 }
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413 void planeelemqt::elem_integration (integratedquant iq, long lcid,long eid,long ri,long ci,vector &nv,vector &x,vector &y)
01414 {
01415 long i,ii,ipp;
01416 double jac,thick;
01417 ivector nodes(nne);
01418 vector w,gp1,gp2,t(nne),ipv,contr(ndofe);
01419 matrix gm;
01420
01421 Mt->give_elemnodes (eid,nodes);
01422 Mc->give_thickness (eid,nodes,t);
01423
01424 fillv (0.0,nv);
01425
01426 for (ii=0;ii<nb;ii++){
01427 allocv (intordsm[ii][ii],gp1);
01428 allocv (intordsm[ii][ii],gp2);
01429 allocv (intordsm[ii][ii],w);
01430 allocm (ncomp[ii],ndofe,gm);
01431 allocv (ncomp[ii],ipv);
01432
01433 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01434 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01435
01436 for (i=0;i<intordsm[ii][ii];i++){
01437 thick = approx (gp1[i],gp2[i],t);
01438
01439 Mm->givequantity (iq,lcid,ipp,cncomp[ii],ipv);
01440
01441
01442 geom_matrix_block (gm,ii,x,y,gp1[i],gp2[i],jac);
01443
01444
01445 mtxv (gm,ipv,contr);
01446
01447 cmulv (jac*w[i]*thick,contr);
01448
01449
01450 addv(contr,nv,nv);
01451
01452 ipp++;
01453 }
01454 destrv (ipv); destrm (gm); destrv (w); destrv (gp1); destrv (gp2);
01455 }
01456 }
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469 void planeelemqt::compute_nlstress (long lcid,long eid,long ri,long ci)
01470 {
01471 long i,ii,ipp;
01472 double **stra;
01473 vector w,gp1,gp2,t(nne),eps(tncomp);
01474
01475 stra = new double* [nne];
01476 for (i=0;i<nne;i++){
01477 stra[i] = new double [tncomp];
01478 }
01479 elem_strains (stra,lcid,eid,ri,ci);
01480
01481 for (ii=0;ii<nb;ii++){
01482 allocv (intordsm[ii][ii],gp1);
01483 allocv (intordsm[ii][ii],gp2);
01484 allocv (intordsm[ii][ii],w);
01485 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01486 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01487
01488 for (i=0;i<intordsm[ii][ii];i++){
01489 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
01490 Mm->storestrain (lcid,ipp,eps);
01491
01492 if (Mp->strcomp==1)
01493 Mm->computenlstresses (ipp);
01494 ipp++;
01495 }
01496 }
01497 for (i=0;i<nne;i++){
01498 delete [] stra[i];
01499 }
01500 delete [] stra;
01501 }
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515 void planeelemqt::local_values (long lcid,long eid,long ri,long ci)
01516 {
01517 long i,ii,ipp;
01518 double **stra;
01519 vector w,gp1,gp2,t(nne),eps(tncomp);
01520
01521 stra = new double* [nne];
01522 for (i=0;i<nne;i++){
01523 stra[i] = new double [tncomp];
01524 }
01525 elem_strains (stra,lcid,eid,ri,ci);
01526
01527 for (ii=0;ii<nb;ii++){
01528 allocv (intordsm[ii][ii],gp1);
01529 allocv (intordsm[ii][ii],gp2);
01530 allocv (intordsm[ii][ii],w);
01531 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01532 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01533
01534 for (i=0;i<intordsm[ii][ii];i++){
01535 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
01536 Mm->storestrain (lcid,ipp,eps);
01537
01538 if (Mp->strcomp==1)
01539 Mm->computenlstresses (ipp);
01540 ipp++;
01541 }
01542 }
01543 for (i=0;i<nne;i++){
01544 delete [] stra[i];
01545 }
01546 delete [] stra;
01547 }
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560 void planeelemqt::compute_nonloc_nlstress (long lcid,long eid,long ri,long ci)
01561 {
01562 long i,ii,ipp;
01563 double **stra;
01564 vector w,gp1,gp2,t(nne),eps(tncomp);
01565
01566 stra = new double* [nne];
01567 for (i=0;i<nne;i++){
01568 stra[i] = new double [tncomp];
01569 }
01570 elem_strains (stra,lcid,eid,ri,ci);
01571
01572 for (ii=0;ii<nb;ii++){
01573 allocv (intordsm[ii][ii],gp1);
01574 allocv (intordsm[ii][ii],gp2);
01575 allocv (intordsm[ii][ii],w);
01576 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01577 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01578
01579 for (i=0;i<intordsm[ii][ii];i++){
01580 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
01581 Mm->storestrain (lcid,ipp,eps);
01582
01583 if (Mp->strcomp==1)
01584 Mm->compnonloc_nlstresses (ipp);
01585 ipp++;
01586 }
01587 }
01588 for (i=0;i<nne;i++){
01589 delete [] stra[i];
01590 }
01591 delete [] stra;
01592 }
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605 void planeelemqt::compute_nlstressincr (long lcid,long eid,long ri,long ci)
01606 {
01607 long i,ii,ipp;
01608 double **stra;
01609 vector w,gp1,gp2,t(nne),eps(tncomp);
01610
01611 stra = new double* [nne];
01612 for (i=0;i<nne;i++){
01613 stra[i] = new double [tncomp];
01614 }
01615 elem_strains (stra,lcid,eid,ri,ci);
01616
01617 for (ii=0;ii<nb;ii++){
01618 allocv (intordsm[ii][ii],gp1);
01619 allocv (intordsm[ii][ii],gp2);
01620 allocv (intordsm[ii][ii],w);
01621 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][ii]);
01622 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01623
01624 for (i=0;i<intordsm[ii][ii];i++){
01625 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
01626 Mm->storestrain (lcid,ipp,eps);
01627
01628 if (Mp->strcomp==1)
01629 Mm->computenlstressesincr (ipp);
01630 ipp++;
01631 }
01632 }
01633 for (i=0;i<nne;i++){
01634 delete [] stra[i];
01635 }
01636 delete [] stra;
01637 }
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650 void planeelemqt::compute_eigstress (long lcid,long eid,long ri,long ci)
01651 {
01652 long i,ii,ipp;
01653 vector eigstr(tncomp),sig(tncomp);
01654 matrix d(tncomp,tncomp);
01655
01656 for (ii=0;ii<nb;ii++){
01657 ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
01658
01659 for (i=0;i<intordsm[ii][ii];i++){
01660 Mm->giveeigstrain (ipp,eigstr);
01661
01662 Mm->matstiff (d,ipp);
01663 mxv (d,eigstr,sig);
01664 Mm->storeeigstress (ipp,sig);
01665 ipp++;
01666 }
01667 }
01668 }
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682 void planeelemqt::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
01683 {
01684 long i,ii;
01685 vector x(nne),y(nne),w(intordsm[ri][ci]),gp1(intordsm[ri][ci]),gp2(intordsm[ri][ci]);
01686
01687 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ri][ci]);
01688 Mt->give_node_coord2d (x,y,eid);
01689 ii=Mt->elements[eid].ipp[ri][ci];
01690
01691 for (i=0;i<intordsm[ri][ci];i++){
01692 if (ii==ipp){
01693 coord[0]=approx (gp1[i],gp2[i],x);
01694 coord[1]=approx (gp1[i],gp2[i],y);
01695 coord[2]=0.0;
01696 }
01697 ii++;
01698 }
01699 }
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711 void planeelemqt::ipcoordblock (long eid,long ri,long ci,double **coord)
01712 {
01713 vector x(nne),y(nne),w(intordsm[ri][ci]),gp1(intordsm[ri][ci]),gp2(intordsm[ri][ci]);
01714
01715 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ri][ci]);
01716 Mt->give_node_coord2d (x,y,eid);
01717
01718 for (long i=0;i<intordsm[ri][ci];i++){
01719 coord[i][0]=approx (gp1[i],gp2[i],x);
01720 coord[i][1]=approx (gp1[i],gp2[i],y);
01721 coord[i][2]=0.0;
01722 }
01723 }
01724
01725 void planeelemqt::nodeforces (long eid,long *le,double *nv,vector &nf)
01726 {
01727 long i;
01728 double xi,eta,jac;
01729 vector x(nne),y(nne),gp(intordb),w(intordb),av(ndofe),v(ndofe);
01730 matrix n(napfun,ndofe),am(ndofe,ndofe);
01731
01732 Mt->give_node_coord2d (x,y,eid);
01733 gauss_points (gp.a,w.a,intordb);
01734
01735 if (le[0]==1){
01736 fillm (0.0,am);
01737 for (i=0;i<intordb;i++){
01738
01739 xi=(1.0-gp[i])/2.0; eta=(1.0+gp[i])/2.0;
01740
01741 bf_matrix (n,xi,eta);
01742
01743 jac1d_2d (jac,x,y,gp[i],0);
01744
01745 jac*=w[i];
01746
01747 nnj (am.a,n.a,jac,n.m,n.n);
01748 }
01749 fillv (0.0,av);
01750 av[0]=nv[0]; av[1]=nv[1]; av[6]=nv[2]; av[7]=nv[3]; av[2]=nv[4]; av[3]=nv[5];
01751 mxv (am,av,v); addv (nf,v,nf);
01752 }
01753 if (le[1]==1){
01754 fillm (0.0,am);
01755 for (i=0;i<intordb;i++){
01756
01757 xi=0.0; eta=(1.0-gp[i])/2.0;
01758
01759 bf_matrix (n,xi,eta);
01760
01761 jac1d_2d (jac,x,y,gp[i],1);
01762
01763 jac*=w[i];
01764
01765 nnj (am.a,n.a,jac,n.m,n.n);
01766 }
01767 fillv (0.0,av);
01768 av[2]=nv[6]; av[3]=nv[7]; av[8]=nv[8]; av[9]=nv[9]; av[4]=nv[10]; av[5]=nv[11];
01769 mxv (am,av,v); addv (nf,v,nf);
01770 }
01771 if (le[2]==1){
01772 fillm (0.0,am);
01773 for (i=0;i<intordb;i++){
01774 xi=(1.0+gp[i])/2.0; eta=0.0;
01775
01776 bf_matrix (n,xi,eta);
01777
01778 jac1d_2d (jac,x,y,gp[i],2);
01779
01780 jac*=w[i];
01781
01782 nnj (am.a,n.a,jac,n.m,n.n);
01783 }
01784 fillv (0.0,av);
01785 av[4]=nv[12]; av[5]=nv[13]; av[10]=nv[14]; av[11]=nv[15]; av[0]=nv[16]; av[1]=nv[17];
01786 mxv (am,av,v); addv (nf,v,nf);
01787 }
01788 }
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814 void planeelemqt::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
01815 {
01816 long i, j, k, ipp;
01817 long ii, jj, nv = nodval.n;
01818 double xi, eta, ipval;
01819 vector w, gp1, gp2, anv(nne);
01820 long nstra, nstre, ncompstr, ncompeqother;
01821 long idstra, idstre, idoth, idic;
01822 inictype ict;
01823
01824 nstra = idstra = nstre = idstre = idoth = idic = 0;
01825
01826 ict = ictn[0];
01827 for (i=0; i<nne; i++)
01828 {
01829 if (ictn[i] != ict)
01830 {
01831 print_err("Incompatible types of initial conditions on element %ld\n"
01832 " at %ld. and %ld. nodes", __FILE__, __LINE__, __func__, eid+1, 1, i+1);
01833 abort();
01834 }
01835 }
01836 for (j = 0; j < nv; j++)
01837 {
01838 for(i = 0; i < nne; i++)
01839 anv[i] = nodval[i][j];
01840 for (ii = 0; ii < nb; ii++)
01841 {
01842 for (jj = 0; jj < nb; jj++)
01843 {
01844 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01845 if (intordsm[ii][jj] == 0)
01846 continue;
01847 allocv (intordsm[ii][jj],gp1);
01848 allocv (intordsm[ii][jj],gp2);
01849 allocv (intordsm[ii][jj],w);
01850 gauss_points_tr (gp1.a, gp2.a, w.a, intordsm[ii][jj]);
01851 for (k = 0; k < intordsm[ii][jj]; k++)
01852 {
01853 xi=gp1[k];
01854 eta=gp2[k];
01855
01856 ipval = approx (xi,eta,anv);
01857 ncompstr = Mm->ip[ipp].ncompstr;
01858 ncompeqother = Mm->ip[ipp].ncompeqother;
01859 if ((ictn[0] & inistrain) && (j < ncompstr))
01860 {
01861 Mm->ip[ipp].strain[idstra] += ipval;
01862 ipp++;
01863 continue;
01864 }
01865 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01866 {
01867 Mm->ip[ipp].stress[idstre] += ipval;
01868 ipp++;
01869 continue;
01870 }
01871 if ((ictn[0] & iniother) && (j < nstra+nstre+ncompeqother))
01872 {
01873 Mm->ip[ipp].eqother[idoth] += ipval;
01874 ipp++;
01875 continue;
01876 }
01877 if ((ictn[0] & inicond) && (j < nv))
01878 {
01879 if (Mm->ic[ipp] == NULL)
01880 {
01881 Mm->ic[ipp] = new double[nv-j];
01882 memset(Mm->ic[ipp], 0, sizeof(*Mm->ic[ipp])*(nv-j));
01883 }
01884 Mm->ic[ipp][idic] += ipval;
01885 ipp++;
01886 continue;
01887 }
01888 ipp++;
01889 }
01890 destrv(gp1); destrv (gp2); destrv (w);
01891 }
01892 }
01893 ipp=Mt->elements[eid].ipp[ri][ci];
01894 ncompstr = Mm->ip[ipp].ncompstr;
01895 ncompeqother = Mm->ip[ipp].ncompeqother;
01896 if ((ictn[0] & inistrain) && (j < ncompstr))
01897 {
01898 nstra++;
01899 idstra++;
01900 continue;
01901 }
01902 if ((ictn[0] & inistress) && (j < nstra + ncompstr))
01903 {
01904 nstre++;
01905 idstre++;
01906 continue;
01907 }
01908 if ((ictn[0] & iniother) && (j < nstra + nstre + ncompeqother))
01909 {
01910 idoth++;
01911 continue;
01912 }
01913 if ((ictn[0] & inicond) && (j < nv))
01914 {
01915 idic++;
01916 continue;
01917 }
01918 }
01919 }
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929 void planeelemqt::ipvolume (long eid,long ri,long ci)
01930 {
01931 long i,ii,jj,ipp;
01932 double jac,thick;
01933 ivector nodes(nne);
01934 vector x(nne),y(nne),t(nne),gp1,gp2,w;
01935
01936
01937 Mt->give_elemnodes (eid,nodes);
01938
01939 Mc->give_thickness (eid,nodes,t);
01940 Mt->give_node_coord2d (x,y,eid);
01941
01942 for (ii=0;ii<nb;ii++){
01943 for (jj=0;jj<nb;jj++){
01944 if (intordsm[ii][jj]==0) continue;
01945
01946 allocv (intordsm[ii][jj],gp1);
01947 allocv (intordsm[ii][jj],gp2);
01948 allocv (intordsm[ii][jj],w);
01949
01950 gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[ii][jj]);
01951 ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
01952
01953 for (i=0;i<intordsm[ii][jj];i++){
01954
01955 thick = approx (gp1[i],gp2[i],t);
01956 jac_2d (jac,x,y,gp1[i],gp2[i]);
01957 jac*=w[i]*thick;
01958
01959 Mm->storeipvol (ipp,jac);
01960
01961 ipp++;
01962 }
01963 destrv (gp1); destrv (gp2); destrv (w);
01964 }
01965 }
01966 }
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982 void planeelemqt::ntdbr_vector (long eid,vector &ntdbr)
01983 {
01984 long intord = 4;
01985 long i,k,l,ipp,ri,ci,lcid;
01986 double thick,jac,**stra;
01987 ivector nodes(nne);
01988 vector x(nne),y(nne),gp1(intord),gp2(intord),w(intord);
01989 vector t(nne),sig(tncomp),eps(tncomp),bf(nne);
01990 matrix d(tncomp,tncomp);
01991
01992 lcid = ri = ci = 0;
01993 ipp = Mt->elements[eid].ipp[ri][ci];
01994
01995 stra = new double* [nne];
01996 for (i=0;i<nne;i++){
01997 stra[i] = new double [tncomp];
01998 }
01999 elem_strains (stra,lcid,eid,ri,ci);
02000
02001 Mt->give_elemnodes (eid,nodes);
02002 Mt->give_node_coord2d (x,y,eid);
02003 Mc->give_thickness (eid,nodes,t);
02004 Mm->matstiff (d,ipp);
02005
02006 gauss_points_tr (gp1.a,gp2.a,w.a,intord);
02007
02008 fillv(0.0,ntdbr);
02009
02010 for (i=0;i<intord;i++){
02011 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
02012 mxv (d,eps,sig);
02013
02014 bf_quad_3_2d (bf.a,gp1[i],gp2[i]);
02015
02016 thick = approx (gp1[i],gp2[i],t);
02017 jac_2d (jac,x,y,gp1[i],gp2[i]);
02018 jac *= w[i]*thick;
02019
02020 for (k=0;k<tncomp;k++)
02021 for (l=0;l<nne;l++)
02022 ntdbr[k*nne+l] += jac * bf[l] * sig[k];
02023
02024 }
02025
02026 for (i=0;i<nne;i++) delete [] stra[i];
02027 delete [] stra;
02028 }
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040 void planeelemqt::ntn_matrix (long eid,matrix &ntn)
02041 {
02042 long intord = 6;
02043 long i,k,l;
02044 double thick,jac;
02045 ivector nodes(nne);
02046 vector x(nne),y(nne),t(nne),gp1(intord),gp2(intord),w(intord),bf(nne);
02047
02048 Mt->give_elemnodes (eid,nodes);
02049 Mt->give_node_coord2d (x,y,eid);
02050 Mc->give_thickness (eid,nodes,t);
02051
02052 gauss_points_tr (gp1.a,gp2.a,w.a,intord);
02053
02054 fillm(0.0,ntn);
02055
02056 for (i=0;i<intord;i++){
02057 bf_quad_3_2d (bf.a,gp1[i],gp2[i]);
02058
02059 thick = approx (gp1[i],gp2[i],t);
02060 jac_2d (jac,x,y,gp1[i],gp2[i]);
02061 jac *= w[i]*thick;
02062
02063 for (k=0;k<nne;k++)
02064 for (l=0;l<nne;l++)
02065 ntn[k][l] += jac * bf[k] * bf[l];
02066
02067 }
02068 }
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093 double planeelemqt :: compute_error (long eid, double &e2, double &u2, double &sizel, vector *rsigfull)
02094 {
02095 long intorde = 7;
02096 long i,ipp,ri,ci,lcid;
02097 double area,thick,jac,contr,**stra;
02098 double zero=1.0e-20;
02099 ivector nodes(nne);
02100 vector x(nne),y(nne),t(nne),gp1(intorde),gp2(intorde),w(intorde),bf(nne);
02101 vector sig_star(tncomp),sig_roof(tncomp),sig_err(tncomp),eps(tncomp);
02102 matrix d(tncomp,tncomp),dinv(tncomp,tncomp);
02103
02104 ri = ci = lcid = 0;
02105 ipp = Mt->elements[eid].ipp[ri][ci];
02106
02107 stra = new double* [nne];
02108 for (i=0;i<nne;i++){
02109 stra[i] = new double [tncomp];
02110 }
02111 elem_strains (stra,lcid,eid,ri,ci);
02112
02113 Mt->give_elemnodes (eid,nodes);
02114 Mt->give_node_coord2d (x,y,eid);
02115 Mc->give_thickness (eid,nodes,t);
02116
02117 Mm->matstiff (d,ipp);
02118 invm (d,dinv,zero);
02119
02120 gauss_points_tr (gp1.a,gp2.a,w.a,intorde);
02121
02122 e2 = u2 = 0;
02123 for (i=0;i<intorde;i++){
02124 thick = approx (gp1[i],gp2[i],t);
02125 jac_2d (jac,x,y,gp1[i],gp2[i]);
02126 jac *= w[i]*thick;
02127
02128 bf_quad_3_2d (bf.a,gp1[i],gp2[i]);
02129 give_der_star (bf,rsigfull,nodes,sig_star,Mt->nn);
02130
02131 appval (gp1[i],gp2[i],0,tncomp,eps,stra);
02132 mxv (d,eps,sig_roof);
02133
02134 subv (sig_star,sig_roof,sig_err);
02135
02136 vxmxv (sig_err,dinv,contr);
02137 e2 += jac*contr;
02138
02139 vxmxv (eps,d,contr);
02140 u2 += jac*contr;
02141 }
02142
02143 area = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
02144 sizel = sqrt(1.1547005384*area);
02145
02146 for (i=0;i<nne;i++) delete [] stra[i];
02147 delete [] stra;
02148
02149 return area/2.0;
02150 }
02151
02152
02153
02154
02155
02156
02157 void planeelemqt::give_sig_roof (matrix &d,double *areacoord,double *x,double *y,long *etnodes,double **xyrr,vector &sig_roof,vector &sig_star)
02158 {
02159 long i,nnecm;
02160 double xx,yy,detcm,shear,xicm,etacm,auxd;
02161 vector areacoordcm(3),eps(3),xcm,ycm,rcm,nodsig_x,nodsig_y,nodsig_xy;
02162 matrix gm;
02163
02164 nnecm = etnodes[1];
02165
02166 allocv (nnecm,xcm);
02167 allocv (nnecm,ycm);
02168 allocv (nnecm,nodsig_x);
02169 allocv (nnecm,nodsig_y);
02170 allocv (nnecm,nodsig_xy);
02171 allocv (2*nnecm,rcm);
02172 allocm (3,2*nnecm,gm);
02173
02174 for (i=0;i<nnecm;i++){
02175 xcm[i] = xyrr[etnodes[i+2]][0];
02176 ycm[i] = xyrr[etnodes[i+2]][1];
02177 rcm[2*i] = xyrr[etnodes[i+2]][2];
02178 rcm[2*i+1] = xyrr[etnodes[i+2]][3];
02179 nodsig_x[i] = xyrr[etnodes[i+2]][4];
02180 nodsig_y[i] = xyrr[etnodes[i+2]][5];
02181 nodsig_xy[i] = xyrr[etnodes[i+2]][6];
02182 }
02183
02184 xx = yy = 0.0;
02185 for (i=0;i<3;i++){
02186 xx += areacoord[i] * x[i];
02187 yy += areacoord[i] * y[i];
02188 }
02189
02190 switch (etnodes[0]){
02191 case planeelementlt:{
02192 detcm = (xcm[1]-xcm[0])*(ycm[2]-ycm[0])-(xcm[2]-xcm[0])*(ycm[1]-ycm[0]);
02193 areacoordcm[0] = ( (xcm[1]*ycm[2] - xcm[2]*ycm[1]) + (ycm[1] - ycm[2])*xx + (xcm[2] - xcm[1])*yy)/detcm;
02194 areacoordcm[1] = ( (xcm[2]*ycm[0] - xcm[0]*ycm[2]) + (ycm[2] - ycm[0])*xx + (xcm[0] - xcm[2])*yy)/detcm;
02195 areacoordcm[2] = ( (xcm[0]*ycm[1] - xcm[1]*ycm[0]) + (ycm[0] - ycm[1])*xx + (xcm[1] - xcm[0])*yy)/detcm;
02196
02197 Pelt->geom_matrix (gm,xcm,ycm);
02198 mxv (gm,rcm,eps);
02199 mxv (d,eps,sig_roof);
02200
02201 sig_star[0] = Pelt->approx (areacoordcm,nodsig_x);
02202 sig_star[1] = Pelt->approx (areacoordcm,nodsig_y);
02203 sig_star[2] = Pelt->approx (areacoordcm,nodsig_xy);
02204 break;
02205 }
02206 case planeelementlq:{
02207 nc_lin_4_2d (1e-6,xx,yy,xcm.a,ycm.a,xicm,etacm);
02208
02209 Pelq->geom_matrix (gm,xcm,ycm,0.0,0.0,auxd);
02210 mxv (gm,rcm,eps);
02211 mxv (d,eps,sig_roof);
02212 shear=sig_roof[2];
02213
02214 Pelq->geom_matrix (gm,xcm,ycm,xicm,etacm,auxd);
02215 mxv (gm,rcm,eps);
02216 mxv (d,eps,sig_roof);
02217 sig_roof[2]=shear;
02218
02219 sig_star[0] = Pelq->approx (xicm,etacm,nodsig_x);
02220 sig_star[1] = Pelq->approx (xicm,etacm,nodsig_y);
02221 sig_star[2] = Pelq->approx (xicm,etacm,nodsig_xy);
02222 break;
02223 }
02224 case planeelementqq:{
02225 nc_lin_4_2d (1e-6,xx,yy,xcm.a,ycm.a,xicm,etacm);
02226
02227 Peqq->geom_matrix (gm,xcm,ycm,xicm,etacm,auxd);
02228 mxv (gm,rcm,eps);
02229 mxv (d,eps,sig_roof);
02230
02231 sig_star[0] = Peqq->approx (xicm,etacm,nodsig_x);
02232 sig_star[1] = Peqq->approx (xicm,etacm,nodsig_y);
02233 sig_star[2] = Peqq->approx (xicm,etacm,nodsig_xy);
02234 break;
02235 }
02236 case planeelementqt:{
02237 detcm = (xcm[1]-xcm[0])*(ycm[2]-ycm[0])-(xcm[2]-xcm[0])*(ycm[1]-ycm[0]);
02238 areacoordcm[0] = ( (xcm[1]*ycm[2] - xcm[2]*ycm[1]) + (ycm[1] - ycm[2])*xx + (xcm[2] - xcm[1])*yy)/detcm;
02239 areacoordcm[1] = ( (xcm[2]*ycm[0] - xcm[0]*ycm[2]) + (ycm[2] - ycm[0])*xx + (xcm[0] - xcm[2])*yy)/detcm;
02240 areacoordcm[2] = ( (xcm[0]*ycm[1] - xcm[1]*ycm[0]) + (ycm[0] - ycm[1])*xx + (xcm[1] - xcm[0])*yy)/detcm;
02241
02242
02243 mxv (gm,rcm,eps);
02244 mxv (d,eps,sig_roof);
02245
02246 sig_star[0] = approx (xicm,etacm,nodsig_x);
02247 sig_star[1] = approx (xicm,etacm,nodsig_y);
02248 sig_star[2] = approx (xicm,etacm,nodsig_xy);
02249 break;
02250 }
02251 default:{
02252 fprintf (stderr,"\n\n unknown element type is required in function");
02253 fprintf (stderr," planeelemqt::give_sig_roof (%s, line %d)",__FILE__,__LINE__);
02254 }
02255 }
02256 }
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322 void planeelemqt :: elchar (long eid, matrix &spsig)
02323 {
02324 long intord = 3;
02325 long i,ipp,ri,ci;
02326 vector eps(tncomp);
02327 matrix d(tncomp,tncomp);
02328
02329 allocm(intord, tncomp, spsig);
02330
02331 ri = ci = 0;
02332 ipp = Mt->elements[eid].ipp[ri][ci];
02333 Mm->matstiff (d,ipp);
02334
02335
02336
02337
02338 for (i=0; i<intord; i++) {
02339 eps[0] = Mm->ip[ipp ].strain[0];
02340 eps[1] = Mm->ip[ipp ].strain[1];
02341 eps[2] = Mm->ip[ipp+3].strain[2];
02342 ipp++;
02343 mxv (d.a, eps.a, spsig.a + i*tncomp, d.m, d.n);
02344 }
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368 }
02369