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