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