00001
00002
00003
00004
00005
00006
00007 #include "globalt.h"
00008 #include "quadlineart.h"
00009 #include "genfile.h"
00010 #include "globmatt.h"
00011
00012 quadlineart::quadlineart (void)
00013 {
00014 long i;
00015
00016
00017 nne=4;
00018
00019 ned=4;
00020
00021 nned=2;
00022
00023 ncomp=2;
00024
00025
00026 ntm=Tp->ntm;
00027
00028
00029 dofe = new long* [ntm];
00030 nip = new long* [ntm];
00031 intordkm = new long* [ntm];
00032 intordcm = new long* [ntm];
00033 ordering = new long* [ntm];
00034 for (i=0;i<ntm;i++){
00035 dofe[i] = new long [ntm];
00036 nip[i] = new long [ntm];
00037 intordkm[i] = new long [ntm];
00038 intordcm[i] = new long [ntm];
00039 ordering[i] = new long [nne];
00040 }
00041
00042
00043 switch (Tp->tmatt){
00044 case nomedium:{ break; }
00045 case onemedium:{
00046 ordering[0][0]=1; ordering[0][1]=2; ordering[0][2]=3; ordering[0][3]=4;
00047 dofe[0][0]=4; intordkm[0][0]=2; intordcm[0][0]=2;
00048 nip[0][0]=8;
00049 ndofe=4; napfun=1;
00050 break;
00051 }
00052 case twomediacoup:{
00053 ordering[0][0]=1; ordering[0][1]=3; ordering[0][2]=5; ordering[0][3]=7;
00054 ordering[1][0]=2; ordering[1][1]=4; ordering[1][2]=6; ordering[1][3]=8;
00055
00056 intordkm[0][0]=2; intordkm[0][1]=2; intordkm[1][0]=2; intordkm[1][1]=2;
00057 intordcm[0][0]=2; intordcm[0][1]=2; intordcm[1][0]=2; intordcm[1][1]=2;
00058
00059 if (Tp->savemode==0){
00060 nip[0][0]=8; nip[0][1]=8; nip[1][0]=8; nip[1][1]=8;
00061 }
00062 if (Tp->savemode==1){
00063 nip[0][0]=8; nip[0][1]=0; nip[1][0]=0; nip[1][1]=0;
00064 }
00065
00066
00067 dofe[0][0]=4; dofe[0][1]=4; dofe[1][0]=4; dofe[1][1]=4;
00068 ndofe=8; napfun=2;
00069 break;
00070 }
00071 case threemediacoup:{
00072 ordering[0][0]=1; ordering[0][1] =4; ordering[0][2] =7; ordering[0][3] =10;
00073 ordering[1][0]=2; ordering[1][1] =5; ordering[1][2] =8; ordering[1][3] =11;
00074 ordering[2][0]=3; ordering[2][1] =6; ordering[2][2] =9; ordering[2][3] =12;
00075
00076 intordkm[0][0]=2; intordkm[0][1]=2; intordkm[0][2]=2;
00077 intordkm[1][0]=2; intordkm[1][1]=2; intordkm[1][2]=2;
00078 intordkm[2][0]=2; intordkm[2][1]=2; intordkm[2][2]=2;
00079
00080 intordcm[0][0]=2; intordcm[0][1]=2; intordcm[0][2]=2;
00081 intordcm[1][0]=2; intordcm[1][1]=2; intordcm[1][2]=2;
00082 intordcm[2][0]=2; intordcm[2][1]=2; intordcm[2][2]=2;
00083
00084 if (Tp->savemode==0){
00085 nip[0][0]=8; nip[0][1]=8; nip[0][2]=8;
00086 nip[1][0]=8; nip[1][1]=8; nip[1][2]=8;
00087 nip[2][0]=8; nip[2][1]=8; nip[2][2]=8;
00088 }
00089 if (Tp->savemode==1){
00090 nip[0][0]=8; nip[0][1]=0; nip[0][2]=0;
00091 nip[1][0]=0; nip[1][1]=0; nip[1][2]=0;
00092 nip[2][0]=0; nip[2][1]=0; nip[2][2]=0;
00093 }
00094
00095 dofe[0][0]=4; dofe[0][1]=4; dofe[0][2]=4;
00096 dofe[1][0]=4; dofe[1][1]=4; dofe[1][2]=4;
00097 dofe[2][0]=4; dofe[2][1]=4; dofe[2][2]=4;
00098
00099
00100 ndofe=12; napfun=3;
00101 break;
00102 }
00103 case fourmediacoup:{
00104 ordering[0][0]=1; ordering[0][1] =5; ordering[0][2] =9; ordering[0][3] =13;
00105 ordering[1][0]=2; ordering[1][1] =6; ordering[1][2] =10; ordering[1][3] =14;
00106 ordering[2][0]=3; ordering[2][1] =7; ordering[2][2] =11; ordering[2][3] =15;
00107 ordering[3][0]=4; ordering[3][1] =8; ordering[3][2] =12; ordering[3][3] =16;
00108
00109 intordkm[0][0]=2; intordkm[0][1]=2; intordkm[0][2]=2; intordkm[0][3]=2;
00110 intordkm[1][0]=2; intordkm[1][1]=2; intordkm[1][2]=2; intordkm[1][3]=2;
00111 intordkm[2][0]=2; intordkm[2][1]=2; intordkm[2][2]=2; intordkm[2][3]=2;
00112 intordkm[3][0]=2; intordkm[3][1]=2; intordkm[3][2]=2; intordkm[3][3]=2;
00113
00114 intordcm[0][0]=2; intordcm[0][1]=2; intordcm[0][2]=2; intordcm[0][3]=2;
00115 intordcm[1][0]=2; intordcm[1][1]=2; intordcm[1][2]=2; intordcm[1][3]=2;
00116 intordcm[2][0]=2; intordcm[2][1]=2; intordcm[2][2]=2; intordcm[2][3]=2;
00117 intordcm[3][0]=2; intordcm[3][1]=2; intordcm[3][2]=2; intordcm[3][3]=2;
00118
00119 if (Tp->savemode==0){
00120 nip[0][0]=8; nip[0][1]=8; nip[0][2]=8; nip[0][3]=8;
00121 nip[1][0]=8; nip[1][1]=8; nip[1][2]=8; nip[1][3]=8;
00122 nip[2][0]=8; nip[2][1]=8; nip[2][2]=8; nip[2][3]=8;
00123 nip[3][0]=8; nip[3][1]=8; nip[3][2]=8; nip[3][3]=8;
00124 }
00125 if (Tp->savemode==1){
00126 nip[0][0]=8; nip[0][1]=0; nip[0][2]=0; nip[0][3]=0;
00127 nip[1][0]=0; nip[1][1]=0; nip[1][2]=0; nip[1][3]=0;
00128 nip[2][0]=0; nip[2][1]=0; nip[2][2]=0; nip[2][3]=0;
00129 nip[3][0]=0; nip[3][1]=0; nip[3][2]=0; nip[3][3]=0;
00130 }
00131
00132 dofe[0][0]=4; dofe[0][1]=4; dofe[0][2]=4; dofe[0][3]=4;
00133 dofe[1][0]=4; dofe[1][1]=4; dofe[1][2]=4; dofe[1][3]=4;
00134 dofe[2][0]=4; dofe[2][1]=4; dofe[2][2]=4; dofe[2][3]=4;
00135 dofe[3][0]=4; dofe[3][1]=4; dofe[3][2]=4; dofe[3][3]=4;
00136
00137 ndofe=16; napfun=4;
00138 break;
00139 }
00140 default:{
00141 print_err("unknown number of transported matters is required",__FILE__,__LINE__,__func__);
00142 }
00143 }
00144 }
00145
00146 quadlineart::~quadlineart (void)
00147 {
00148 long i;
00149
00150 for (i=0;i<ntm;i++){
00151 delete [] dofe[i];
00152 delete [] nip[i];
00153 delete [] intordkm[i];
00154 delete [] intordcm[i];
00155 delete [] ordering[i];
00156 }
00157 delete [] dofe;
00158 delete [] nip;
00159 delete [] intordkm;
00160 delete [] intordcm;
00161 delete [] ordering;
00162 }
00163
00164
00165
00166
00167
00168
00169
00170
00171 void quadlineart::codnum (long *cn,long ri)
00172 {
00173 long i;
00174 for (i=0;i<nne;i++){
00175 cn[i]=ordering[ri][i];
00176 }
00177 }
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 double quadlineart::approx (double xi,double eta,vector &nodval)
00188 {
00189 double f;
00190 vector bf(nne);
00191
00192 bf_lin_4_2d (bf.a,xi,eta);
00193
00194 scprd (bf,nodval,f);
00195
00196 return f;
00197 }
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211 void quadlineart::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
00212 {
00213 long i,j,ii;
00214 double xi,eta;
00215 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]);
00216
00217 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00218 Tt->give_node_coord2d (x,y,eid);
00219
00220 if (Tp->savemode==0)
00221 ii=Tt->elements[eid].ipp[ri][ci];
00222 if (Tp->savemode==1)
00223 ii=Tt->elements[eid].ipp[0][0];
00224
00225
00226 for (i=0;i<intordkm[ri][ci];i++){
00227 xi=gp[i];
00228 for (j=0;j<intordkm[ri][ci];j++){
00229 eta=gp[j];
00230 if (ii==ipp){
00231 coord[0]=approx (xi,eta,x);
00232 coord[1]=approx (xi,eta,y);
00233 coord[2]=0.0;
00234 }
00235 ii++;
00236 }
00237 }
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 void quadlineart::ipcoordblock (long eid,long ri,long ci,double **coord)
00251 {
00252 long i,j,k;
00253 double xi,eta;
00254 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]);
00255
00256 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00257 Tt->give_node_coord2d (x,y,eid);
00258
00259 k=0;
00260 for (i=0;i<intordkm[ri][ci];i++){
00261 xi=gp[i];
00262 for (j=0;j<intordkm[ri][ci];j++){
00263 eta=gp[j];
00264
00265 coord[k][0]=approx (xi,eta,x);
00266 coord[k][1]=approx (xi,eta,y);
00267 coord[k++][2]=0.0;
00268 }
00269 }
00270 }
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 void quadlineart::intpointval (long eid)
00281 {
00282 long i,j,k,ii,jj,ipp;
00283 double xi,eta,val;
00284 ivector cn(ndofe);
00285 vector r(ndofe),t(nne),gp,w;
00286
00287 Tt->give_code_numbers (eid,cn.a);
00288
00289 elemvalues (0,eid,r.a,cn.a,ndofe);
00290
00291 for (k=0;k<Tp->ntm;k++){
00292
00293
00294 for (i=0;i<dofe[k][k];i++){
00295 t[i]=r[ordering[k][i]-1];
00296 }
00297
00298 for (ii=0;ii<Tp->ntm;ii++){
00299 for (jj=0;jj<Tp->ntm;jj++){
00300
00301
00302 reallocv (intordkm[ii][jj],gp);
00303 reallocv (intordkm[ii][jj],w);
00304 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00305
00306 ipp=Tt->elements[eid].ipp[ii][jj];
00307 for (i=0;i<intordkm[ii][jj];i++){
00308 xi=gp[i];
00309 for (j=0;j<intordkm[ii][jj];j++){
00310 eta=gp[j];
00311
00312 val = approx (xi,eta,t);
00313 Tm->ip[ipp].av[k]=val;
00314 ipp++;
00315 }
00316 }
00317
00318
00319 reallocv (intordcm[ii][jj],gp);
00320 reallocv (intordcm[ii][jj],w);
00321 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00322
00323 for (i=0;i<intordcm[ii][jj];i++){
00324 xi=gp[i];
00325 for (j=0;j<intordcm[ii][jj];j++){
00326 eta=gp[j];
00327
00328 val = approx (xi,eta,t);
00329 Tm->ip[ipp].av[k]=val;
00330 ipp++;
00331 }
00332 }
00333
00334 if (Tp->savemode==1) break;
00335 }
00336 if (Tp->savemode==1) break;
00337 }
00338 }
00339
00340 }
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 void quadlineart::intpointval (long eid,vector &nodval,vector &ipval)
00357 {
00358 long i,j,kk,ii,jj,ipp;
00359 double xi,eta;
00360 vector gp,w;
00361
00362 kk = 0;
00363
00364 for (ii=0;ii<Tp->ntm;ii++){
00365 for (jj=0;jj<Tp->ntm;jj++){
00366
00367
00368 reallocv (intordkm[ii][jj],gp);
00369 reallocv (intordkm[ii][jj],w);
00370 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00371
00372 ipp=Tt->elements[eid].ipp[ii][jj];
00373 for (i=0;i<intordkm[ii][jj];i++){
00374 xi=gp[i];
00375 for (j=0;j<intordkm[ii][jj];j++){
00376 eta=gp[j];
00377
00378 ipval[kk] = approx (xi,eta,nodval);
00379
00380 ipp++;
00381 kk++;
00382 }
00383 }
00384
00385
00386
00387 reallocv (intordcm[ii][jj],gp);
00388 reallocv (intordcm[ii][jj],w);
00389 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00390
00391 for (i=0;i<intordcm[ii][jj];i++){
00392 xi=gp[i];
00393 for (j=0;j<intordcm[ii][jj];j++){
00394 eta=gp[j];
00395
00396 ipval[kk] = approx (xi,eta,nodval);
00397
00398 ipp++;
00399 }
00400 }
00401
00402 if (Tp->savemode==1) break;
00403 }
00404 if (Tp->savemode==1) break;
00405 }
00406 }
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417 void quadlineart::intpointval_puc (long eid)
00418 {
00419 long i,j,k,ii,jj,ipp;
00420 double xi,eta,val;
00421 ivector cn(ndofe);
00422 vector r(ndofe),t(nne),gp,w;
00423
00424 Tt->give_code_numbers (eid,cn.a);
00425
00426 elemvalues_puc (0,eid,r.a,ndofe);
00427
00428 for (k=0;k<Tp->ntm;k++){
00429
00430
00431 for (i=0;i<dofe[k][k];i++){
00432 t[i]=r[ordering[k][i]-1];
00433 }
00434
00435 for (ii=0;ii<Tp->ntm;ii++){
00436 for (jj=0;jj<Tp->ntm;jj++){
00437
00438
00439 reallocv (intordkm[ii][jj],gp);
00440 reallocv (intordkm[ii][jj],w);
00441 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00442
00443 ipp=Tt->elements[eid].ipp[ii][jj];
00444 for (i=0;i<intordkm[ii][jj];i++){
00445 xi=gp[i];
00446 for (j=0;j<intordkm[ii][jj];j++){
00447 eta=gp[j];
00448
00449 val = approx (xi,eta,t);
00450 Tm->ip[ipp].av[k]=val;
00451 ipp++;
00452 }
00453 }
00454
00455
00456
00457 reallocv (intordcm[ii][jj],gp);
00458 reallocv (intordcm[ii][jj],w);
00459 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00460
00461 for (i=0;i<intordcm[ii][jj];i++){
00462 xi=gp[i];
00463 for (j=0;j<intordcm[ii][jj];j++){
00464 eta=gp[j];
00465
00466 val = approx (xi,eta,t);
00467 Tm->ip[ipp].av[k]=val;
00468 ipp++;
00469 }
00470 }
00471
00472 if (Tp->savemode==1) break;
00473 }
00474 if (Tp->savemode==1) break;
00475 }
00476 }
00477
00478 }
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488 void quadlineart::intpointgrad (long eid)
00489 {
00490 long i,j,ii,jj,k,ipp;
00491 double xi,eta,jac;
00492 ivector cn(ndofe);
00493 vector x(nne),y(nne),r(ndofe),t(nne),gp,w,grad(ncomp);
00494 matrix gm(ncomp,nne);
00495
00496
00497 Tt->give_node_coord2d (x,y,eid);
00498
00499 Tt->give_code_numbers (eid,cn.a);
00500
00501 elemvalues (0,eid,r.a,cn.a,ndofe);
00502
00503 for (k=0;k<Tp->ntm;k++){
00504
00505
00506 for (i=0;i<dofe[k][k];i++){
00507 t[i]=r[ordering[k][i]-1];
00508 }
00509
00510 for (ii=0;ii<Tp->ntm;ii++){
00511 for (jj=0;jj<Tp->ntm;jj++){
00512
00513
00514 reallocv (intordkm[ii][jj],gp);
00515 reallocv (intordkm[ii][jj],w);
00516 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00517
00518 ipp=Tt->elements[eid].ipp[ii][jj];
00519 for (i=0;i<intordkm[ii][jj];i++){
00520 xi=gp[i];
00521 for (j=0;j<intordkm[ii][jj];j++){
00522 eta=gp[j];
00523
00524
00525 grad_matrix (gm,x,y,xi,eta,jac);
00526 mxv (gm,t,grad);
00527 Tm->storegrad (k,ipp,grad);
00528 ipp++;
00529 }
00530 }
00531
00532
00533 reallocv (intordcm[ii][jj],gp);
00534 reallocv (intordcm[ii][jj],w);
00535 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00536
00537 for (i=0;i<intordcm[ii][jj];i++){
00538 xi=gp[i];
00539 for (j=0;j<intordcm[ii][jj];j++){
00540 eta=gp[j];
00541
00542
00543 grad_matrix (gm,x,y,xi,eta,jac);
00544 mxv (gm,t,grad);
00545 Tm->storegrad (k,ipp,grad);
00546 ipp++;
00547 }
00548 }
00549
00550 if (Tp->savemode==1) break;
00551 }
00552 if (Tp->savemode==1) break;
00553 }
00554 }
00555
00556 }
00557
00558
00559
00560
00561
00562
00563
00564
00565 void quadlineart::intpointother (long eid)
00566 {
00567 long i,j,k,ii,jj,ipp,ncompo,nodid;
00568 double xi,eta,val,*r;
00569 ivector nodes(nne);
00570 vector t(nne),gp,w;
00571
00572
00573 Tt->give_elemnodes (eid,nodes);
00574
00575
00576 nodid=nodes[0];
00577
00578
00579 ncompo=Tt->nodes[nodid].ncompother;
00580 r = new double [ncompo*nne];
00581
00582
00583 nodalotherval (r,nodes);
00584
00585 for (k=0;k<ncompo;k++){
00586
00587
00588 for (i=0;i<nne;i++){
00589 t[i]=r[i*ncompo+k];
00590 }
00591
00592 for (ii=0;ii<Tp->ntm;ii++){
00593 for (jj=0;jj<Tp->ntm;jj++){
00594
00595
00596 reallocv (intordkm[ii][jj],gp);
00597 reallocv (intordkm[ii][jj],w);
00598 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00599
00600 ipp=Tt->elements[eid].ipp[ii][jj];
00601 for (i=0;i<intordkm[ii][jj];i++){
00602 xi=gp[i];
00603 for (j=0;j<intordkm[ii][jj];j++){
00604 eta=gp[j];
00605
00606 val = approx (xi,eta,t);
00607 Tm->ip[ipp].other[k]=val;
00608 ipp++;
00609 }
00610 }
00611
00612
00613
00614 reallocv (intordcm[ii][jj],gp);
00615 reallocv (intordcm[ii][jj],w);
00616 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00617
00618 for (i=0;i<intordcm[ii][jj];i++){
00619 xi=gp[i];
00620 for (j=0;j<intordcm[ii][jj];j++){
00621 eta=gp[j];
00622
00623 val = approx (xi,eta,t);
00624 Tm->ip[ipp].other[k]=val;
00625 ipp++;
00626 }
00627 }
00628
00629 if (Tp->savemode==1) break;
00630 }
00631 if (Tp->savemode==1) break;
00632 }
00633 }
00634
00635 delete [] r;
00636 }
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646 void quadlineart::bf_matrix (matrix &n,double xi,double eta)
00647 {
00648 fillm (0.0,n);
00649 bf_lin_4_2d (n.a,xi,eta);
00650 }
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662 void quadlineart::grad_matrix (matrix &gm,vector &x,vector &y,double xi,double eta,double &jac)
00663 {
00664 long i;
00665 vector dx(nne),dy(nne);
00666
00667 dx_bf_lin_4_2d (dx.a,eta);
00668 dy_bf_lin_4_2d (dy.a,xi);
00669
00670 derivatives_2d (dx,dy,jac,x,y,xi,eta);
00671
00672 fillm (0.0,gm);
00673
00674 for (i=0;i<nne;i++){
00675 gm[0][i]=dx[i];
00676 gm[1][i]=dy[i];
00677 }
00678
00679 }
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694 void quadlineart::conductivity_matrix (long lcid,long eid,long ri,long ci,matrix &km)
00695 {
00696 long i,j,ii;
00697 double thick,xi,eta,ww1,ww2,jac;
00698 ivector nodes(nne);
00699 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
00700 matrix gm(ncomp,dofe[ri][ci]),d;
00701
00702 matrix n(1,dofe[ri][ci]);
00703
00704
00705 Tt->give_elemnodes (eid,nodes);
00706
00707 Tc->give_thickness (eid,nodes,t);
00708
00709 Tt->give_node_coord2d (x,y,eid);
00710
00711 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00712
00713 fillm (0.0,km);
00714
00715 if (Tp->savemode==0)
00716 ii=Tt->elements[eid].ipp[ri][ci];
00717 if (Tp->savemode==1)
00718 ii=Tt->elements[eid].ipp[0][0];
00719
00720 for (i=0;i<intordkm[ri][ci];i++){
00721 xi=gp[i]; ww1=w[i];
00722 for (j=0;j<intordkm[ri][ci];j++){
00723 eta=gp[j]; ww2=w[j];
00724
00725
00726 grad_matrix (gm,x,y,xi,eta,jac);
00727
00728
00729 reallocm(ncomp,ncomp,d);
00730 Tm->matcond (d,ii,ri,ci);
00731
00732
00733 thick = approx (xi,eta,t);
00734
00735 if (jac<0.0){
00736 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00737 jac=fabs(jac);
00738 }
00739
00740 jac*=thick*ww1*ww2;
00741
00742
00743 bdbj (km.a,gm.a,d.a,jac,gm.m,gm.n);
00744
00745
00746 reallocm(1,ncomp,d);
00747 Tm->matcond2(d,ii,ri,ci);
00748 bf_matrix (n, xi, eta);
00749 bdbjac(km, n, d, gm, jac);
00750
00751 ii++;
00752 }
00753 }
00754
00755
00756 if ((Tt->elements[eid].transi[lcid]==3)||(Tt->elements[eid].transi[lcid]==4)){
00757
00758 transmission_matrix (lcid,eid,ri,ci,km);
00759 }
00760
00761 }
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776 void quadlineart::l_matrix (long lcid,long eid,long ri,long ci,matrix &lm)
00777 {
00778 long i,j,ii;
00779 double thick,xi,eta,ww1,ww2,jac;
00780 ivector nodes(nne);
00781 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
00782 matrix gm(ncomp,dofe[ri][ci]),d;
00783
00784
00785 Tt->give_elemnodes (eid,nodes);
00786
00787 Tc->give_thickness (eid,nodes,t);
00788
00789 Tt->give_node_coord2d (x,y,eid);
00790
00791 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00792
00793 fillm (0.0,lm);
00794
00795 if (Tp->savemode==0)
00796 ii=Tt->elements[eid].ipp[ri][ci];
00797 if (Tp->savemode==1)
00798 ii=Tt->elements[eid].ipp[0][0];
00799
00800 for (i=0;i<intordkm[ri][ci];i++){
00801 xi=gp[i]; ww1=w[i];
00802 for (j=0;j<intordkm[ri][ci];j++){
00803 eta=gp[j]; ww2=w[j];
00804
00805
00806 grad_matrix (gm,x,y,xi,eta,jac);
00807
00808
00809 reallocm(ncomp,ncomp,d);
00810 Tm->matcond (d,ii,ri,ci);
00811
00812
00813 thick = approx (xi,eta,t);
00814
00815 if (jac<0.0){
00816 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00817 jac=fabs(jac);
00818 }
00819
00820 jac*=thick*ww1*ww2;
00821
00822
00823 mxm(d,gm,lm);
00824 cmulm(jac,lm);
00825
00826 ii++;
00827 }
00828 }
00829 }
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 void quadlineart::l_t_matrix (long lcid,long eid,long ri,long ci,matrix &lm)
00845 {
00846 long i,j,ii;
00847 double thick,xi,eta,ww1,ww2,jac;
00848 ivector nodes(nne);
00849 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
00850 matrix gm(ncomp,dofe[ri][ci]),d;
00851
00852
00853 Tt->give_elemnodes (eid,nodes);
00854
00855 Tc->give_thickness (eid,nodes,t);
00856
00857 Tt->give_node_coord2d (x,y,eid);
00858
00859 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00860
00861 fillm (0.0,lm);
00862
00863 if (Tp->savemode==0)
00864 ii=Tt->elements[eid].ipp[ri][ci];
00865 if (Tp->savemode==1)
00866 ii=Tt->elements[eid].ipp[0][0];
00867
00868 for (i=0;i<intordkm[ri][ci];i++){
00869 xi=gp[i]; ww1=w[i];
00870 for (j=0;j<intordkm[ri][ci];j++){
00871 eta=gp[j]; ww2=w[j];
00872
00873
00874 grad_matrix (gm,x,y,xi,eta,jac);
00875
00876
00877 reallocm(ncomp,ncomp,d);
00878 Tm->matcond (d,ii,ri,ci);
00879
00880
00881 thick = approx (xi,eta,t);
00882
00883 if (jac<0.0){
00884 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00885 jac=fabs(jac);
00886 }
00887
00888 jac*=thick*ww1*ww2;
00889
00890
00891 mtxm(gm,d,lm);
00892 cmulm(jac,lm);
00893
00894 ii++;
00895 }
00896 }
00897 }
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911 void quadlineart::capacity_matrix (long eid,long ri,long ci,matrix &cm)
00912 {
00913 long i,j,ii;
00914 double jac,xi,eta,w1,w2,thick,rho,c;
00915 ivector nodes(nne);
00916 vector x(nne),y(nne),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),t(nne),dens(nne);
00917 matrix n(1,dofe[ri][ci]);
00918
00919 Tt->give_elemnodes (eid,nodes);
00920 Tc->give_thickness (eid,nodes,t);
00921 Tc->give_density (eid,nodes,dens);
00922
00923 Tt->give_node_coord2d (x,y,eid);
00924 gauss_points (gp.a,w.a,intordcm[ri][ci]);
00925
00926 fillm (0.0,cm);
00927
00928 if (Tp->savemode==0)
00929 ii=Tt->elements[eid].ipp[ri][ci]+intordkm[ri][ci]*intordkm[ri][ci];
00930 if (Tp->savemode==1)
00931 ii=Tt->elements[eid].ipp[0][0]+intordkm[0][0]*intordkm[0][0];
00932
00933 for (i=0;i<intordcm[ri][ci];i++){
00934 xi=gp[i]; w1=w[i];
00935 for (j=0;j<intordcm[ri][ci];j++){
00936 eta=gp[j]; w2=w[j];
00937 jac_2d (jac,x,y,xi,eta);
00938 bf_matrix (n,xi,eta);
00939
00940 c=Tm->capcoeff (ii,ri,ci);
00941 thick = approx (xi,eta,t);
00942 rho = approx (xi,eta,dens);
00943 jac*=w1*w2*thick*rho*c;
00944
00945 if (jac<0.0){
00946 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00947 jac=fabs(jac);
00948 }
00949
00950 nnj (cm.a,n.a,jac,n.m,n.n);
00951 ii++;
00952 }
00953 }
00954
00955 }
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970 void quadlineart::quantity_source_vector (vector &sv,vector &nodval,long eid,long ri,long ci)
00971 {
00972 long i,j;
00973 double jac,xi,eta,w1,w2,thick;
00974 ivector nodes(nne);
00975 vector x(nne),y(nne),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),t(nne),dens(nne),v(dofe[ri][ci]);
00976 matrix n(1,dofe[ri][ci]),nm(dofe[ri][ci],dofe[ri][ci]);
00977
00978 Tt->give_elemnodes (eid,nodes);
00979
00980 Tc->give_thickness (eid,nodes,t);
00981
00982 Tt->give_node_coord2d (x,y,eid);
00983 gauss_points (gp.a,w.a,intordcm[ri][ci]);
00984
00985 fillm (0.0,nm);
00986
00987 for (i=0;i<intordcm[ri][ci];i++){
00988 xi=gp[i]; w1=w[i];
00989 for (j=0;j<intordcm[ri][ci];j++){
00990 eta=gp[j]; w2=w[j];
00991 jac_2d (jac,x,y,xi,eta);
00992
00993 bf_matrix (n,xi,eta);
00994
00995 thick = approx (xi,eta,t);
00996
00997 if (jac<0.0){
00998 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
00999 jac=fabs(jac);
01000 }
01001
01002 jac*=w1*w2*thick;
01003
01004 nnj (nm.a,n.a,jac,n.m,n.n);
01005 }
01006 }
01007 mxv (nm,nodval,v); addv (sv,v,sv);
01008
01009 }
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022 void quadlineart::internal_fluxes (long lcid,long eid,vector &ifl)
01023 {
01024 long i,j,ipp;
01025 double xi,eta,jac,thick;
01026 ivector nodes(nne);
01027 vector x(nne),y(nne),w,gp,t(nne),fl(ncomp),contr(dofe[lcid][lcid]);
01028 matrix gm(ncomp,dofe[lcid][lcid]);
01029
01030
01031 Tt->give_elemnodes (eid,nodes);
01032 Tc->give_thickness (eid,nodes,t);
01033 Tt->give_node_coord2d (x,y,eid);
01034
01035 reallocv (intordkm[lcid][lcid],gp);
01036 reallocv (intordkm[lcid][lcid],w);
01037
01038 gauss_points (gp.a,w.a,intordkm[lcid][lcid]);
01039
01040 if (Tp->savemode==0)
01041 ipp=Tt->elements[eid].ipp[lcid][lcid];
01042 if (Tp->savemode==1)
01043 ipp=Tt->elements[eid].ipp[0][0];
01044
01045
01046 for (i=0;i<intordkm[lcid][lcid];i++){
01047 xi=gp[i];
01048 for (j=0;j<intordkm[lcid][lcid];j++){
01049 eta=gp[j];
01050
01051
01052 thick = approx (xi,eta,t);
01053
01054 Tm->computenlfluxes (lcid,ipp);
01055
01056 Tm->givefluxes (lcid,ipp,fl);
01057
01058 grad_matrix (gm,x,y,xi,eta,jac);
01059
01060 mtxv (gm,fl,contr);
01061
01062 if (jac<0.0){
01063 print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
01064 jac=fabs(jac);
01065 }
01066
01067 cmulv (jac*w[i]*w[j]*thick,contr);
01068
01069 addv (contr,ifl,ifl);
01070
01071 ipp++;
01072 }
01073 }
01074
01075 }
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087 void quadlineart::res_conductivity_matrix (long eid,long lcid,matrix &km)
01088 {
01089 long i,j,*rcn,*ccn;
01090 matrix lkm;
01091
01092 for (i=0;i<ntm;i++){
01093 for (j=0;j<ntm;j++){
01094 rcn = new long [dofe[i][j]];
01095 ccn = new long [dofe[i][j]];
01096 reallocm (dofe[i][j],dofe[i][j],lkm);
01097
01098 conductivity_matrix (i,eid,i,j,lkm);
01099
01100 codnum (rcn,i);
01101 codnum (ccn,j);
01102 mat_localize (km,lkm,rcn,ccn);
01103 delete [] rcn; delete [] ccn;
01104 }
01105 }
01106
01107 }
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121 void quadlineart::res_l_matrix (long eid,long lcid,matrix &lm)
01122 {
01123 long i,j,*rcn,*ccn;
01124 matrix lkm;
01125
01126 for (i=0;i<ntm;i++){
01127 for (j=0;j<ntm;j++){
01128 rcn = new long [2];
01129 ccn = new long [dofe[i][j]];
01130 if (i==0){
01131 rcn[0]=1;
01132 rcn[1]=2;
01133 }
01134 if (i==1){
01135 rcn[0]=3;
01136 rcn[1]=4;
01137 }
01138 reallocm (ncomp,dofe[i][j],lkm);
01139 l_matrix (i,eid,i,j,lkm);
01140 codnum (ccn,j);
01141
01142 mat_localize (lm,lkm,rcn,ccn);
01143 delete [] rcn; delete [] ccn;
01144 }
01145 }
01146 }
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161 void quadlineart::res_l_t_matrix (long eid,long lcid,matrix &lm)
01162 {
01163 long i,j,*rcn,*ccn;
01164 matrix lkm;
01165
01166 for (i=0;i<ntm;i++){
01167 for (j=0;j<ntm;j++){
01168 rcn = new long [dofe[i][j]];
01169 ccn = new long [2];
01170 if (j==0){
01171 ccn[0]=1;
01172 ccn[1]=2;
01173 }
01174 if (j==1){
01175 ccn[0]=3;
01176 ccn[1]=4;
01177 }
01178 reallocm (dofe[i][j],ncomp,lkm);
01179 l_t_matrix (i,eid,i,j,lkm);
01180 codnum (rcn,i);
01181
01182 mat_localize (lm,lkm,rcn,ccn);
01183 delete [] rcn; delete [] ccn;
01184 }
01185 }
01186 }
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199 void quadlineart::averd_matrix (long eid,matrix &lm)
01200 {
01201 long i,j,ii,*rcn,*ccn;
01202 matrix d(ncomp,ncomp);
01203
01204 rcn = new long [ncomp];
01205 ccn = new long [ncomp];
01206
01207 for (i=0;i<ntm;i++){
01208 for (j=0;j<ntm;j++){
01209
01210 if (Tp->savemode==0)
01211 ii=Tt->elements[eid].ipp[i][j];
01212 if (Tp->savemode==1)
01213 ii=Tt->elements[eid].ipp[0][0];
01214
01215
01216 Tm->matcond (d,ii,i,j);
01217
01218 if (i==0){
01219 rcn[0]=1;
01220 rcn[1]=2;
01221 }
01222 if (i==1){
01223 rcn[0]=3;
01224 rcn[1]=4;
01225 }
01226
01227 if (j==0){
01228 ccn[0]=1;
01229 ccn[1]=2;
01230 }
01231 if (j==1){
01232 ccn[0]=3;
01233 ccn[1]=4;
01234 }
01235
01236 mat_localize (lm,d,rcn,ccn);
01237 }
01238 }
01239
01240 delete [] rcn; delete [] ccn;
01241 }
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255 void quadlineart::averc_matrix (long eid,matrix &lm)
01256 {
01257 long i,j,ii;
01258 double c;
01259
01260 for (i=0;i<ntm;i++){
01261 for (j=0;j<ntm;j++){
01262
01263 if (Tp->savemode==0)
01264 ii=Tt->elements[eid].ipp[i][j];
01265 if (Tp->savemode==1)
01266 ii=Tt->elements[eid].ipp[0][0];
01267
01268
01269 c = Tm->capcoeff (ii,i,j);
01270
01271 lm[i][j] = c*Tm->ip[ii].av[j];
01272 }
01273 }
01274 }
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284 double quadlineart::elem_area (long eid)
01285 {
01286 double area,jac;
01287 vector x(nne),y(nne);
01288 matrix gm(ncomp,dofe[0][0]);
01289
01290 Tt->give_node_coord2d (x,y,eid);
01291
01292
01293 grad_matrix (gm,x,y,0.0,1.0,jac);
01294
01295
01296 area = jac*4.0;
01297
01298 return area;
01299 }
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312 void quadlineart::res_capacity_matrix (long eid,matrix &cm)
01313 {
01314 long i,j,*rcn,*ccn,k,l;
01315 double s;
01316 matrix lcm;
01317
01318 for (i=0;i<ntm;i++){
01319 for (j=0;j<ntm;j++){
01320 rcn = new long [dofe[i][j]];
01321 ccn = new long [dofe[i][j]];
01322 reallocm (dofe[i][j],dofe[i][j],lcm);
01323
01324 capacity_matrix (eid,i,j,lcm);
01325
01326
01327 if (Tp->diagcap == 1){
01328 for (k=0;k<dofe[i][j];k++){
01329 s=0.0;
01330 for (l=0;l<dofe[i][j];l++){
01331 s+=lcm[k][l];
01332 lcm[k][l]=0.0;
01333 }
01334 lcm[k][k]=s;
01335 }
01336 }
01337
01338 codnum (rcn,i);
01339 codnum (ccn,j);
01340
01341 mat_localize (cm,lcm,rcn,ccn);
01342 delete [] rcn; delete [] ccn;
01343 }
01344 }
01345 }
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357 void quadlineart::res_convection_vector (vector &f,long lcid,long eid,long leid)
01358 {
01359 long *cn;
01360 vector lf;
01361
01362
01363
01364 if ((Tt->elements[eid].transi[lcid]==2)||(Tt->elements[eid].transi[lcid]==4)){
01365
01366 cn = new long [dofe[lcid][lcid]];
01367
01368 reallocv (dofe[lcid][lcid],lf);
01369
01370 convection_vector (lf,lcid,eid,leid,lcid);
01371
01372
01373 codnum (cn,lcid);
01374
01375 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
01376
01377 delete [] cn;
01378
01379 cmulv (-1.0,f,f);
01380 }
01381 }
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393 void quadlineart::res_transmission_vector (vector &f,long lcid,long eid,long leid)
01394 {
01395 long *cn,i;
01396 vector lf;
01397
01398 cn = new long [dofe[lcid][lcid]];
01399
01400 codnum (cn,lcid);
01401 reallocv (dofe[lcid][lcid],lf);
01402
01403
01404
01405 if ((Tt->elements[eid].transi[lcid]==3)||(Tt->elements[eid].transi[lcid]==4)){
01406
01407 for (i=0;i<ntm;i++){
01408 nullv (lf.a,lf.n);
01409
01410 transmission_vector (lf,lcid,eid,leid,i);
01411
01412 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
01413 }
01414 }
01415 delete [] cn;
01416 }
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428 void quadlineart::res_quantity_source_vector (vector &sv,vector &nodval,long lcid,long eid)
01429 {
01430 long *cn;
01431 vector lsv;
01432
01433 cn = new long [dofe[lcid][lcid]];
01434 reallocv (dofe[lcid][lcid],lsv);
01435 quantity_source_vector (lsv,nodval,eid,lcid,lcid);
01436 codnum (cn,lcid);
01437 locglob (sv.a,lsv.a,cn,dofe[lcid][lcid]);
01438 delete [] cn;
01439 }
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449 void quadlineart::res_internal_fluxes (long eid,vector &elemif)
01450 {
01451 long i,*cn;
01452 vector lif,tdnv;
01453 matrix cm;
01454
01455 for (i=0;i<ntm;i++){
01456 cn = new long [dofe[i][i]];
01457 reallocv (dofe[i][i],lif);
01458 internal_fluxes (i,eid,lif);
01459 codnum (cn,i);
01460 locglob (elemif.a,lif.a,cn,dofe[i][i]);
01461
01462 delete [] cn;
01463 }
01464
01465 cn = new long [ndofe];
01466 Gtt->give_code_numbers (eid,cn);
01467 reallocv (ndofe,tdnv);
01468 reallocv (ndofe,lif);
01469 reallocm (ndofe,ndofe,cm);
01470 nodalderivatives (tdnv.a,cn,ndofe);
01471 res_capacity_matrix (eid,cm);
01472 mxv (cm,tdnv,lif);
01473
01474 subv (elemif,lif,elemif);
01475
01476 delete [] cn;
01477 }
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489 double quadlineart::total_integral(long eid,vector &nodval)
01490 {
01491 long i,j;
01492 double thick,xi,eta,ww1,ww2,jac,value,f;
01493 ivector nodes(nne);
01494 vector x(nne),y(nne),w(2),gp(2),t(nne);
01495
01496 Tt->give_elemnodes (eid,nodes);
01497 Tc->give_thickness (eid,nodes,t);
01498
01499 Tt->give_node_coord2d (x,y,eid);
01500 gauss_points (gp.a,w.a,2);
01501
01502 f = 0.0;
01503
01504 for (i=0;i<2;i++){
01505 xi=gp[i]; ww1=w[i];
01506 for (j=0;j<2;j++){
01507 eta=gp[j]; ww2=w[j];
01508
01509 jac_2d (jac,x,y,xi,eta);
01510 thick = approx (xi,eta,t);
01511 value = approx (xi,eta,nodval);
01512
01513 jac*=thick*ww1*ww2*value;
01514 f = f + jac;
01515 }
01516 }
01517 return(f);
01518 }
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531 void quadlineart::res_boundary_flux (vector &f,long lcid,long eid,long leid)
01532 {
01533 long *cn,i;
01534 vector lf;
01535
01536 cn = new long [dofe[lcid][lcid]];
01537 codnum (cn,lcid);
01538
01539
01540 for (i=0;i<ntm;i++){
01541 reallocv (dofe[lcid][lcid],lf);
01542
01543 if ((Tt->elements[eid].transi[i]==3)||(Tt->elements[eid].transi[i]==4)){
01544 boundary_flux (lf,i,eid,leid,lcid,i);
01545 }
01546
01547 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
01548 }
01549
01550 delete [] cn;
01551 }
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566 void quadlineart::volume_rhs_vector (long lcid,long eid,long ri,long ci,vector &vrhs)
01567 {
01568 long i,j,ii;
01569 double thick,xi,eta,ww1,ww2,jac;
01570 ivector nodes(nne);
01571 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
01572 matrix gm(ncomp,dofe[ri][ci]),d;
01573 matrix km(dofe[ri][ci],2);
01574
01575
01576 Tt->give_elemnodes (eid,nodes);
01577
01578 Tt->give_node_coord2d (x,y,eid);
01579
01580 gauss_points (gp.a,w.a,intordkm[ri][ci]);
01581
01582 fillm (0.0,km);
01583
01584 if (Tp->savemode==0)
01585 ii=Tt->elements[eid].ipp[ri][ci];
01586 if (Tp->savemode==1)
01587 ii=Tt->elements[eid].ipp[0][0];
01588
01589 for (i=0;i<intordkm[ri][ci];i++){
01590 xi=gp[i]; ww1=w[i];
01591 for (j=0;j<intordkm[ri][ci];j++){
01592 eta=gp[j]; ww2=w[j];
01593
01594 grad_matrix (gm,x,y,xi,eta,jac);
01595
01596
01597 reallocm(ncomp,ncomp,d);
01598
01599 Tm->volume_rhs (d,ii,ri,ci,ncomp);
01600
01601
01602 thick = approx (xi,eta,t);
01603
01604 jac*=ww1*ww2*thick;
01605
01606
01607 nnjac (km,gm,d,jac);
01608
01609 ii++;
01610 }
01611 }
01612 for (i=0;i<vrhs.n;i++){
01613 vrhs[i] = km[i][0];
01614 }
01615 }
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628 void quadlineart::res_volume_rhs_vector (vector &f,long eid,long lcid)
01629 {
01630 long i,*cn;
01631 vector lf;
01632
01633 for (i=0;i<ntm;i++){
01634 cn = new long [dofe[i][i]];
01635 reallocv (dofe[i][i],lf);
01636 codnum (cn,i);
01637 volume_rhs_vector (i,eid,i,i,lf);
01638 locglob (f.a,lf.a,cn,dofe[i][i]);
01639 delete [] cn;
01640 }
01641 }
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658 void quadlineart::convection_vector (vector &v,long lcid,long eid,long leid,long cid)
01659 {
01660 long i,ipp;
01661 bocontypet *bc;
01662 ivector nodes(nne);
01663 vector x(nne),y(nne),w(intordkm[lcid][cid]),gp(intordkm[lcid][cid]),t(nne);
01664 vector list(nned*ned),trc(nned*ned),nodval(nne),av(nne),coef(nne);
01665 matrix km(nne,nne);
01666
01667 fillv (0.0,v);
01668
01669
01670 Tt->give_elemnodes (eid,nodes);
01671
01672 Tc->give_thickness (eid,nodes,t);
01673
01674 Tt->give_node_coord2d (x,y,eid);
01675
01676
01677 gauss_points (gp.a,w.a,intordkm[lcid][cid]);
01678
01679
01680 bc = new bocontypet [ned];
01681 Tb->lc[lcid].elemload[leid].give_bc (bc);
01682
01683 Tb->lc[lcid].elemload[leid].give_nodval (lcid,list);
01684
01685
01686 for (i=0;i<nned*ned;i++){
01687 trc[i]=1.0;
01688 }
01689
01690
01691 if (Tp->savemode==0)
01692 ipp=Tt->elements[eid].ipp[lcid][cid];
01693 if (Tp->savemode==1)
01694 ipp=Tt->elements[eid].ipp[0][0];
01695
01696
01697 for (i=0;i<ned;i++){
01698
01699 if (bc[i]==2 || bc[i]==3 || bc[i]==5){
01700
01701 edgenodeval (i,nodval,list);
01702
01703 edgenodeval (i,coef,trc);
01704
01705 fillm (0.0,km);
01706
01707
01708 edge_integral (i,x,y,intordkm[lcid][cid],gp,w,t,coef,km);
01709
01710 mxv (km,nodval,av);
01711
01712 addv (v,av,v);
01713 }
01714 }
01715
01716 delete [] bc;
01717 }
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731 void quadlineart::transmission_matrix (long lcid,long eid,long ri,long ci,matrix &km)
01732 {
01733 long i,j,ipp,leid;
01734 bocontypet *bc;
01735 ivector nodes(nne);
01736 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]);
01737 vector trc(nned*ned),coeff(nne),t(nne);
01738
01739
01740 Tt->give_elemnodes (eid,nodes);
01741
01742 Tc->give_thickness (eid,nodes,t);
01743
01744 Tt->give_node_coord2d (x,y,eid);
01745
01746
01747 gauss_points (gp.a,w.a,intordkm[ri][ci]);
01748
01749
01750 for (i=0;i<Tb->lc[lcid].neb;i++){
01751 if (Tb->lc[lcid].elemload[i].eid==eid){
01752 leid=i;
01753
01754
01755 bc = new bocontypet [ned];
01756 Tb->lc[lcid].elemload[leid].give_bc (bc);
01757
01758 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
01759
01760
01761 if (Tp->savemode==0)
01762 ipp=Tt->elements[eid].ipp[ri][ci];
01763 if (Tp->savemode==1)
01764 ipp=Tt->elements[eid].ipp[0][0];
01765
01766
01767 for (j=0;j<ned;j++){
01768
01769 if (bc[j]==5 || bc[j]>10){
01770
01771 transf_coeff (j,coeff,trc,eid,ri,ci,ipp,bc,0);
01772
01773 edge_integral (j,x,y,intordkm[ri][ci],gp,w,t,coeff,km);
01774 }
01775 }
01776
01777 }
01778 }
01779
01780 delete [] bc;
01781 }
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797 void quadlineart::transmission_vector (vector &v,long lcid,long eid,long leid,long cid)
01798 {
01799 long i,ipp;
01800 bocontypet *bc;
01801 ivector nodes(nne);
01802 vector x(nne),y(nne),z(nne),w(intordkm[lcid][cid]),gp(intordkm[lcid][cid]),t(nne);
01803 vector list(nned*ned),trc(nned*ned),trr(nned*ned),nodval(nne),av(nne),coef(nne);
01804 matrix km(nne,nne);
01805
01806 fillv (0.0,v);
01807
01808
01809 Tt->give_elemnodes (eid,nodes);
01810
01811 Tc->give_thickness (eid,nodes,t);
01812
01813 Tt->give_node_coord2d (x,y,eid);
01814
01815
01816 gauss_points (gp.a,w.a,intordkm[lcid][cid]);
01817
01818
01819 bc = new bocontypet [ned];
01820 Tb->lc[lcid].elemload[leid].give_bc (bc);
01821
01822 Tb->lc[lcid].elemload[leid].give_external_nodval (lcid,cid,list);
01823
01824 Tb->lc[lcid].elemload[leid].give_trc (lcid,cid,trc);
01825
01826 Tb->lc[lcid].elemload[leid].give_trr (lcid,trr);
01827
01828
01829 if (Tp->savemode==0)
01830 ipp=Tt->elements[eid].ipp[lcid][cid];
01831 if (Tp->savemode==1)
01832 ipp=Tt->elements[eid].ipp[0][0];
01833
01834
01835
01836 for (i=0;i<ned;i++){
01837
01838 if (bc[i]==5 || bc[i]>10){
01839
01840 transf_val (i,nodval,list,trc,trr,eid,lcid,cid,ipp,bc);
01841
01842 transf_coeff (i,coef,trc,eid,lcid,cid,ipp,bc,1);
01843
01844 fillm (0.0,km);
01845
01846 edge_integral (i,x,y,intordkm[lcid][cid],gp,w,t,coef,km);
01847
01848 mxv (km,nodval,av);
01849
01850 addv (v,av,v);
01851 }
01852 }
01853
01854 delete [] bc;
01855 }
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871 void quadlineart::boundary_flux (vector &v,long lcid,long eid,long leid,long ri,long ci)
01872 {
01873 long i,ipp;
01874 bocontypet *bc;
01875 ivector nodes(nne);
01876 vector x(nne),y(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),t(nne);
01877 vector list(nned*ned),trc(nned*ned),trr(nned*ned),nodval(nne),av(nne),coef(nne);
01878 matrix km(nne,nne);
01879
01880 fillv (0.0,v);
01881
01882
01883 Tt->give_elemnodes (eid,nodes);
01884
01885 Tc->give_thickness (eid,nodes,t);
01886
01887 Tt->give_node_coord2d (x,y,eid);
01888
01889
01890 gauss_points (gp.a,w.a,intordkm[ri][ci]);
01891
01892
01893 bc = new bocontypet [ned];
01894 Tb->lc[lcid].elemload[leid].give_bc (bc);
01895
01896 Tb->lc[lcid].elemload[leid].give_nodval (lcid,list);
01897
01898 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
01899
01900 Tb->lc[lcid].elemload[leid].give_trr (lcid,trr);
01901
01902
01903 if (Tp->savemode==0)
01904 ipp=Tt->elements[eid].ipp[ri][ci];
01905 if (Tp->savemode==1)
01906 ipp=Tt->elements[eid].ipp[0][0];
01907
01908 for (i=0;i<ned;i++){
01909
01910 if (bc[i]==5 || bc[i]>10){
01911
01912 transf_flux (i,nodval,list,trc,trr,eid,ri,ci,ipp,bc);
01913
01914 edgenodeval (i,coef,trc);
01915
01916 fillm (0.0,km);
01917
01918 edge_integral (i,x,y,intordkm[ri][ci],gp,w,t,coef,km);
01919
01920 mxv (km,nodval,av);
01921
01922 addv (v,av,v);
01923 }
01924 }
01925
01926 delete [] bc;
01927 }
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942 void quadlineart::edge_integral (long edg,vector &x,vector &y,long intord,vector &gp,vector &w,
01943 vector &t,vector &coef,matrix &km)
01944 {
01945 long i;
01946 double xi,eta,jac,ipval,thick;
01947 matrix n(1,nne);
01948
01949 if (edg==0){
01950 eta=1.0;
01951 for (i=0;i<intord;i++){
01952 xi=gp[i];
01953
01954 bf_matrix (n,xi,eta);
01955 jac1d_2d (jac,x,y,xi,edg);
01956
01957 ipval=approx (xi,eta,coef);
01958 thick=approx (xi,eta,t);
01959
01960 jac*=w[i]*ipval*thick;
01961 nnj (km.a,n.a,jac,n.m,n.n);
01962 }
01963 }
01964
01965 if (edg==1){
01966 xi=-1.0;
01967 for (i=0;i<intord;i++){
01968 eta=gp[i];
01969
01970 bf_matrix (n,xi,eta);
01971 jac1d_2d (jac,x,y,eta,edg);
01972
01973 ipval=approx (xi,eta,coef);
01974 thick=approx (xi,eta,t);
01975
01976 jac*=w[i]*ipval*thick;
01977 nnj (km.a,n.a,jac,n.m,n.n);
01978 }
01979 }
01980
01981 if (edg==2){
01982 eta=-1.0;
01983 for (i=0;i<intord;i++){
01984 xi=gp[i];
01985
01986 bf_matrix (n,xi,eta);
01987 jac1d_2d (jac,x,y,xi,edg);
01988
01989 ipval=approx (xi,eta,coef);
01990 thick=approx (xi,eta,t);
01991
01992 jac*=w[i]*ipval*thick;
01993 nnj (km.a,n.a,jac,n.m,n.n);
01994 }
01995 }
01996
01997 if (edg==3){
01998 xi=1.0;
01999 for (i=0;i<intord;i++){
02000 eta=gp[i];
02001
02002 bf_matrix (n,xi,eta);
02003 jac1d_2d (jac,x,y,eta,edg);
02004
02005 ipval=approx (xi,eta,coef);
02006 thick=approx (xi,eta,t);
02007
02008 jac*=w[i]*ipval*thick;
02009 nnj (km.a,n.a,jac,n.m,n.n);
02010 }
02011 }
02012
02013 }
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028 void quadlineart::transf_flux (long edg,vector &coeff,vector &list,vector &trc,vector &trr,long eid,long ri,long ci,long ipp,bocontypet *bc)
02029 {
02030 long i,j,k;
02031 double tr,new_flux;
02032 ivector nodes(nne),edgenod(nned);
02033
02034
02035 fillv (0.0,coeff);
02036
02037
02038 Tt->give_elemnodes (eid,nodes);
02039
02040
02041 linquadrilat_edgnod (edgenod.a,edg);
02042
02043
02044 j=edg*nned;
02045
02046 for (i=0;i<nned;i++){
02047 if (bc[edg]==90){
02048 tr = trr[j+i]/trc[j+i];
02049 }
02050
02051
02052 k=edgenod[i];
02053
02054 Tm->transmission_flux(new_flux,list[j+i],tr,ri,ci,nodes[k],bc[edg],ipp);
02055
02056 coeff[k]=new_flux;
02057 }
02058
02059 }
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074 void quadlineart::transf_coeff (long edg,vector &coeff,vector &list,long eid,long ri,long ci,long ipp,bocontypet *bc,int flag)
02075 {
02076 long i,j,k;
02077 double new_coeff;
02078 ivector nodes(nne),edgenod(nned);
02079
02080
02081 fillv (0.0,coeff);
02082
02083
02084 Tt->give_elemnodes (eid,nodes);
02085
02086
02087 linquadrilat_edgnod (edgenod.a,edg);
02088
02089
02090 j=edg*nned;
02091
02092 for (i=0;i<nned;i++){
02093
02094 k=edgenod[i];
02095
02096 Tm->transmission_transcoeff(new_coeff,list[j+i],ri,ci,nodes[k],bc[edg],ipp,flag);
02097
02098 coeff[k]=new_coeff;
02099 }
02100
02101 }
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119 void quadlineart::transf_val (long edg,vector &nodval,vector &list,vector &trc,vector &trr,long eid,long ri,long ci,long ipp,bocontypet *bc)
02120 {
02121 long i,j,k;
02122 double tr,new_nodval;
02123 ivector nodes(nne),edgenod(nned);
02124
02125
02126 fillv (0.0,nodval);
02127
02128
02129 Tt->give_elemnodes (eid,nodes);
02130
02131
02132 linquadrilat_edgnod (edgenod.a,edg);
02133
02134
02135 j=edg*nned;
02136
02137
02138 for (i=0;i<nned;i++){
02139 if (bc[edg]==90){
02140 tr = trr[j+i]/trc[j+i];
02141 }
02142
02143
02144 k=edgenod[i];
02145
02146 Tm->transmission_nodval(new_nodval,list[j+i],tr,ri,ci,nodes[k],bc[edg],ipp);
02147
02148
02149 nodval[k]=new_nodval;
02150 }
02151
02152 }
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163 void quadlineart::edgenodeval (long edg,vector &nodval,vector &list)
02164 {
02165 long i,j;
02166 ivector edgenod(nned);
02167
02168 fillv (0.0,nodval);
02169 linquadrilat_edgnod (edgenod.a,edg);
02170
02171 for (i=0;i<nned;i++){
02172 j=edgenod[i];
02173 nodval[j]=list[edg*nned+i];
02174 }
02175 }
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186 void quadlineart::intpointflux (long eid)
02187 {
02188 long i,ii,jj,k,ipp;
02189
02190 for (k=0;k<Tp->ntm;k++){
02191
02192 for (ii=0;ii<Tp->ntm;ii++){
02193 for (jj=0;jj<Tp->ntm;jj++){
02194
02195 ipp=Tt->elements[eid].ipp[ii][jj];
02196
02197 for (i=0;i<intordkm[ii][jj];i++){
02198
02199 if (Tp->fluxcomp==1)
02200 Tm->computenlfluxes (k,ipp);
02201
02202 ipp++;
02203 }
02204
02205 for (i=0;i<intordcm[ii][jj];i++){
02206
02207 if (Tp->fluxcomp==1)
02208 Tm->computenlfluxes (k,ipp);
02209
02210 ipp++;
02211 }
02212 if (Tp->savemode==1) break;
02213 }
02214 if (Tp->savemode==1) break;
02215 }
02216 }
02217 }
02218
02219
02220
02221
02222
02223
02224
02225
02226
02227
02228
02229 void quadlineart::nod_grads_ip (long eid)
02230 {
02231 long i,j,k,ipp;
02232 ivector ipnum(nne),nod(nne);
02233 vector grad(ncomp);
02234
02235
02236
02237 ipp=Tt->elements[eid].ipp[0][0];
02238 nodip_planelq (ipp,intordkm[0][0],ipnum);
02239
02240
02241 Tt->give_elemnodes (eid,nod);
02242
02243 for (i=0;i<nne;i++){
02244 for (k=0;k<Tp->ntm;k++){
02245
02246 Tm->givegrad (k,ipnum[i],grad);
02247
02248
02249 j=nod[i];
02250 Tt->nodes[j].storegrad (k,ncomp,grad);
02251 }
02252 }
02253
02254 }
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264 void quadlineart::nod_fluxes_ip (long eid)
02265 {
02266 long i,j,k,ipp;
02267 ivector ipnum(nne),nod(nne);
02268 vector flux(ncomp);
02269
02270
02271
02272 ipp=Tt->elements[eid].ipp[0][0];
02273 nodip_planelq (ipp,intordkm[0][0],ipnum);
02274
02275
02276 Tt->give_elemnodes (eid,nod);
02277
02278 intpointflux (eid);
02279
02280 for (i=0;i<nne;i++){
02281 for (k=0;k<Tp->ntm;k++){
02282
02283 Tm->givefluxes (k,ipnum[i],flux);
02284
02285
02286 j=nod[i];
02287 Tt->nodes[j].storeflux (k,flux);
02288 }
02289 }
02290
02291 }
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301 void quadlineart::nod_eqother_ip (long eid)
02302 {
02303 long i,j,k,ipp,ncompo;
02304 ivector ipnum(nne),nod(nne);
02305 vector eqother;
02306
02307
02308
02309 ipp=Tt->elements[eid].ipp[0][0];
02310 nodip_planelq (ipp,intordkm[0][0],ipnum);
02311
02312
02313
02314 Tt->give_elemnodes (eid,nod);
02315
02316 for (i=0;i<nne;i++){
02317 ncompo = Tm->givencompeqother (ipnum[i],0);
02318 reallocv (ncompo,eqother);
02319 Tm->giveeqother (ipnum[i],ncompo,eqother.a);
02320
02321
02322 j=nod[i];
02323 for (k=0; k<ncompo; k++)
02324 Tt->nodes[j].storeeqother (k,eqother(k));
02325
02326 }
02327 }
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339 void quadlineart::nod_others_comp (long lcid,long eid,long ri,long ci)
02340 {
02341 long i,j,k,ipp,ncompother,ii;
02342 vector h,r(ndofe);
02343 ivector nod(nne),cn(ndofe);
02344 double other;
02345
02346 Tt->give_elemnodes (eid,nod);
02347 Tt->give_code_numbers (eid,cn.a);
02348 elemvalues (lcid,eid,r.a,cn.a,ndofe);
02349 ncompother = Tm->givencompother();
02350 reallocv (ntm,h);
02351
02352 ii = 0;
02353 for (i=0;i<nne;i++){
02354 if (Tp->savemode==0)
02355 ipp=Tt->elements[eid].ipp[ri][ci];
02356 if (Tp->savemode==1)
02357 ipp=Tt->elements[eid].ipp[0][0];
02358
02359 for (j=0;j<ntm;j++){
02360 h[j] = r[ii];
02361 ii++;
02362 }
02363
02364 for(k=0;k<ncompother;k++){
02365 other = Tm->givecompother(k,ipp,h.a);
02366 Tt->nodes[nod[i]].storeother (k,other);
02367 }
02368 }
02369
02370 }
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383 void quadlineart::higher_to_lower_level (long eid,long *counter,double *buff)
02384 {
02385 long i,ipp;
02386 vector gr(2);
02387
02388
02389 ipp=Tt->elements[eid].ipp[0][0];
02390
02391
02392 for (i=0;i<ntm;i++){
02393
02394
02395 buff[counter[0]]=Tm->ip[ipp].av[i];
02396 counter[0]++;
02397
02398
02399 Tm->givegrad (i,ipp,gr);
02400 buff[counter[0]]=gr[0];
02401 counter[0]++;
02402 buff[counter[0]]=gr[1];
02403 counter[0]++;
02404
02405 }
02406
02407 }
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433 double quadlineart :: compute_error (long eid, vector *rderfull, int mattid, double &e2, double &u2, double &sizel)
02434 {
02435 long intord = 2;
02436 long i, j, ipp;
02437 double xi, eta, contr, jac, area;
02438 ivector nodes(nne),cn;
02439 vector x(nne), y(nne), t(nne), gp(intord), w(intord), bf(nne), r;
02440 vector der(ncomp), der_fine(ncomp), der_err(ncomp), eps(ncomp);
02441 matrix d(ncomp,ncomp),gm;
02442
02443 Tt->give_elemnodes (eid, nodes);
02444 Tt->give_node_coord2d (x, y, eid);
02445 Tc->give_thickness (eid, nodes, t);
02446
02447
02448 if (Tp->savemode==1) ipp = Tt->elements[eid].ipp[0][0];
02449 else { print_err("save mode is not 1", __FILE__, __LINE__, __func__); exit (1); }
02450
02451 Tm->matcond (d, ipp, mattid, mattid);
02452
02453 gauss_points (gp.a, w.a, intord);
02454
02455 e2 = u2 = 0;
02456 for (i=0; i<intord; i++) {
02457 xi=gp[i];
02458 for (j=0; j<intord; j++) {
02459 eta=gp[j];
02460
02461 Tm->givegrad (mattid, ipp, der);
02462 ipp++;
02463 jac_2d (jac, x, y, xi, eta);
02464
02465 bf_lin_4_2d (bf.a, xi, eta);
02466 give_der_star (bf, rderfull, nodes, der_fine, Tt->nn);
02467 subv (der_fine, der, der_err);
02468
02469 jac *= w[i]*w[j] * approx (xi,eta,t);
02470
02471 vxmxv (der, d, contr); u2 += contr * jac;
02472 vxmxv (der_err, d, contr); e2 += contr * jac;
02473 }
02474 }
02475
02476 area = ( (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]) + (x[2]-x[0])*(y[3]-y[0])-(x[3]-x[0])*(y[2]-y[0]) ) / 2.0;
02477 sizel = sqrt(area);
02478 return area;
02479 }
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491
02492
02493
02494
02495
02496 void quadlineart::transq_nodval (long eid,vector &nodval,nonmechquant nmq)
02497 {
02498 long i,ipid;
02499 ivector ipnum(nne);
02500
02501
02502
02503 ipid=Tt->elements[eid].ipp[0][0];
02504 nodip_planelq (ipid,intordkm[0][0], ipnum);
02505
02506 for (i=0;i<nne;i++)
02507 {
02508
02509 nodval[i] = Tm->givetransq(nmq, ipnum[i]);
02510 }
02511 }
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532 void quadlineart::transq_nodval_comp (long eid,vector &nodval, long ncne, long nq, nonmechquant *qt)
02533 {
02534 long i,j,ipid,ncompo;
02535 ivector ipnum(nne), enod(nne);
02536 intpointst ipb;
02537
02538
02539
02540 ipid=Tt->elements[eid].ipp[0][0];
02541
02542
02543 Tt->give_elemnodes(eid, enod);
02544
02545
02546
02547 nodip_planelq (ipid,intordkm[0][0], ipnum);
02548
02549
02550
02551 ipb.copy(Tm->ip[ipid], ntm);
02552
02553
02554
02555 for (i=0;i<ncne;i++)
02556 {
02557
02558 nodalval (0, Tm->ip[ipid].av, enod[i]);
02559
02560
02561
02562
02563 ncompo = Tm->ip[ipnum[i]].ncompeqother;
02564 if (ncompo){
02565
02566 Tm->storeeqother(ipid, ncompo, Tm->ip[ipnum[i]].eqother);
02567 }
02568 ncompo = Tm->ip[ipnum[i]].ncompother;
02569 if (ncompo){
02570
02571 Tm->storeother(ipid, ncompo, Tm->ip[ipnum[i]].other);
02572 }
02573
02574 Tm->mat_aux_values(ipid);
02575
02576
02577 for (j=0; j<nq; j++)
02578 nodval[j*ncne+i] = Tm->givetransq(qt[j], ipid);
02579 }
02580
02581 Tm->ip[ipid].copy(ipb, ntm);
02582 }