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