00001
00002
00003
00004
00005
00006
00007 #include "globalt.h"
00008 #include "linbartax.h"
00009 #include "genfile.h"
00010 #include "globmatt.h"
00011
00012 linbartax::linbartax (void)
00013 {
00014 long i;
00015
00016
00017 nne=2;
00018
00019 ncomp=1;
00020
00021 nen=2;
00022
00023 ned=2;
00024 nned=1;
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 onemedium:{
00045 ordering[0][0]=1; ordering[0][1]=2;
00046 dofe[0][0]=2; intordkm[0][0]=2; intordcm[0][0]=2; nip[0][0]=4;
00047 ndofe=2; napfun=1;
00048 break;
00049 }
00050 case twomediacoup:{
00051 ordering[0][0]=1; ordering[0][1]=3;
00052 ordering[1][0]=2; ordering[1][1]=4;
00053
00054 intordkm[0][0]=2; intordkm[0][1]=2; intordkm[1][0]=2; intordkm[1][1]=2;
00055 intordcm[0][0]=2; intordcm[0][1]=2; intordcm[1][0]=2; intordcm[1][1]=2;
00056
00057 if (Tp->savemode==0){
00058 nip[0][0]=4; nip[0][1]=4; nip[1][0]=4; nip[1][1]=4;
00059 }
00060 if (Tp->savemode==1){
00061 nip[0][0]=4; nip[0][1]=0; nip[1][0]=0; nip[1][1]=0;
00062 }
00063
00064 dofe[0][0]=2; dofe[0][1]=2; dofe[1][0]=2; dofe[1][1]=2;
00065 ndofe=4; napfun=2;
00066 break;
00067 }
00068 case threemediacoup:{
00069 ordering[0][0]=1; ordering[0][1]=4;
00070 ordering[1][0]=2; ordering[1][1]=5;
00071 ordering[2][0]=3; ordering[2][1]=6;
00072
00073 intordkm[0][0]=2; intordkm[0][1]=2; intordkm[0][2]=2;
00074 intordkm[1][0]=2; intordkm[1][1]=2; intordkm[1][2]=2;
00075 intordkm[2][0]=2; intordkm[2][1]=2; intordkm[2][2]=2;
00076
00077 intordcm[0][0]=2; intordcm[0][1]=2; intordcm[0][2]=2;
00078 intordcm[1][0]=2; intordcm[1][1]=2; intordcm[1][2]=2;
00079 intordcm[2][0]=2; intordcm[2][1]=2; intordcm[2][2]=2;
00080
00081 if (Tp->savemode==0){
00082 nip[0][0]=4; nip[0][1]=4; nip[0][2]=4;
00083 nip[1][0]=4; nip[1][1]=4; nip[1][2]=4;
00084 nip[2][0]=4; nip[2][1]=4; nip[2][2]=4;
00085 }
00086 if (Tp->savemode==1){
00087 nip[0][0]=4; nip[0][1]=0; nip[0][2]=0;
00088 nip[1][0]=0; nip[1][1]=0; nip[1][2]=0;
00089 nip[2][0]=0; nip[2][1]=0; nip[2][2]=0;
00090 }
00091
00092 dofe[0][0]=2; dofe[0][1]=2; dofe[0][2]=2;
00093 dofe[1][0]=2; dofe[1][1]=2; dofe[1][2]=2;
00094 dofe[2][0]=2; dofe[2][1]=2; dofe[2][2]=2;
00095
00096 ndofe=6; napfun=3;
00097 break;
00098 }
00099 default:{
00100 print_err("unknown type of transported matter is required",__FILE__,__LINE__,__func__);
00101 }
00102 }
00103 }
00104
00105 linbartax::~linbartax (void)
00106 {
00107 long i;
00108
00109 for (i=0;i<ntm;i++){
00110 delete [] dofe[i];
00111 delete [] nip[i];
00112 delete [] intordkm[i];
00113 delete [] intordcm[i];
00114 delete [] ordering[i];
00115 }
00116 delete [] dofe;
00117 delete [] nip;
00118 delete [] intordkm;
00119 delete [] intordcm;
00120 delete [] ordering;
00121 }
00122
00123 void linbartax::codnum (long *cn,long ri)
00124 {
00125 long i;
00126 for (i=0;i<nne;i++){
00127 cn[i]=ordering[ri][i];
00128 }
00129 }
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139 double linbartax::approx (double xi,vector &nodval)
00140 {
00141 double f;
00142 vector bf(nne);
00143
00144 bf_lin_1d (bf.a,xi);
00145
00146 scprd (bf,nodval,f);
00147
00148 return f;
00149 }
00150
00151
00152
00153
00154
00155
00156
00157
00158 void linbartax::intpointval (long eid)
00159 {
00160 long i,k,ii,jj,ipp;
00161 double xi,val;
00162 ivector cn(ndofe);
00163 vector r(ndofe),t(nne),gp,w;
00164
00165 Tt->give_code_numbers (eid,cn.a);
00166 elemvalues (0,eid,r.a,cn.a,ndofe);
00167
00168
00169 for (k=0;k<ntm;k++){
00170
00171 for (i=0;i<dofe[k][k];i++){
00172 t[i]=r[ordering[k][i]-1];
00173 }
00174
00175 for (ii=0;ii<Tp->ntm;ii++){
00176 for (jj=0;jj<Tp->ntm;jj++){
00177
00178 reallocv (intordkm[ii][jj],gp);
00179 reallocv (intordkm[ii][jj],w);
00180 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00181
00182 ipp=Tt->elements[eid].ipp[ii][jj];
00183 for (i=0;i<intordkm[ii][jj];i++){
00184 xi=gp[i];
00185
00186 val = approx (xi,t);
00187 Tm->ip[ipp].av[k]=val;
00188 ipp++;
00189 }
00190
00191 reallocv (intordcm[ii][jj],gp);
00192 reallocv (intordcm[ii][jj],w);
00193 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00194
00195 for (i=0;i<intordcm[ii][jj];i++){
00196 xi=gp[i];
00197
00198 val = approx (xi,t);
00199 Tm->ip[ipp].av[k]=val;
00200 ipp++;
00201 }
00202
00203 if (Tp->savemode==1) break;
00204 }
00205 if (Tp->savemode==1) break;
00206 }
00207 }
00208 }
00209
00210
00211
00212
00213
00214
00215
00216
00217 void linbartax::intpointgrad (long eid)
00218 {
00219 long i,ii,jj,k,ipp;
00220 double xi,jac;
00221 ivector cn(ndofe);
00222 vector x(nne),y(nne),r(ndofe),t(nne),gp,w,grad(ncomp);
00223 matrix gm(ncomp,nne);
00224
00225 Tt->give_node_coord2d (x,y,eid);
00226 Tt->give_code_numbers (eid,cn.a);
00227 elemvalues (0,eid,r.a,cn.a,ndofe);
00228
00229 for (k=0;k<Tp->ntm;k++){
00230
00231 for (i=0;i<dofe[k][k];i++){
00232 t[i]=r[ordering[k][i]-1];
00233 }
00234
00235 for (ii=0;ii<Tp->ntm;ii++){
00236 for (jj=0;jj<Tp->ntm;jj++){
00237
00238 reallocv (intordkm[ii][jj],gp);
00239 reallocv (intordkm[ii][jj],w);
00240 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00241
00242 ipp=Tt->elements[eid].ipp[ii][jj];
00243 for (i=0;i<intordkm[ii][jj];i++){
00244 xi=gp[i];
00245
00246
00247 grad_matrix (gm,x,xi,jac);
00248 mxv (gm,t,grad);
00249 Tm->storegrad (k,ipp,grad);
00250 ipp++;
00251 }
00252
00253 reallocv (intordcm[ii][jj],gp);
00254 reallocv (intordcm[ii][jj],w);
00255 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00256
00257 for (i=0;i<intordcm[ii][jj];i++){
00258 xi=gp[i];
00259
00260
00261 grad_matrix (gm,x,xi,jac);
00262 mxv (gm,t,grad);
00263 Tm->storegrad (k,ipp,grad);
00264 ipp++;
00265 }
00266
00267 if (Tp->savemode==1) break;
00268 }
00269 if (Tp->savemode==1) break;
00270 }
00271 }
00272 }
00273
00274
00275
00276
00277
00278
00279
00280
00281 void linbartax::intpointother (long eid)
00282 {
00283 long i,k,ii,jj,ipp,ncompo,nodid;
00284 double xi,val,*r;
00285 ivector nodes(nne);
00286 vector t(nne),gp,w;
00287
00288
00289 Tt->give_elemnodes (eid,nodes);
00290
00291
00292 nodid=nodes[0];
00293
00294
00295 ncompo=Tt->nodes[nodid].ncompother;
00296 r = new double [ncompo*nne];
00297
00298
00299 nodalotherval (r,nodes);
00300
00301 for (k=0;k<ncompo;k++){
00302
00303 for (i=0;i<nne;i++){
00304 t[i]=r[i*ncompo+k];
00305 }
00306
00307 for (ii=0;ii<Tp->ntm;ii++){
00308 for (jj=0;jj<Tp->ntm;jj++){
00309
00310 reallocv (intordkm[ii][jj],gp);
00311 reallocv (intordkm[ii][jj],w);
00312 gauss_points (gp.a,w.a,intordkm[ii][jj]);
00313
00314 ipp=Tt->elements[eid].ipp[ii][jj];
00315 for (i=0;i<intordkm[ii][jj];i++){
00316 xi=gp[i];
00317
00318 val = approx (xi,t);
00319 Tm->ip[ipp].other[k]=val;
00320 ipp++;
00321 }
00322
00323 reallocv (intordcm[ii][jj],gp);
00324 reallocv (intordcm[ii][jj],w);
00325 gauss_points (gp.a,w.a,intordcm[ii][jj]);
00326
00327 for (i=0;i<intordcm[ii][jj];i++){
00328 xi=gp[i];
00329
00330 val = approx (xi,t);
00331 Tm->ip[ipp].other[k]=val;
00332 ipp++;
00333 }
00334
00335 if (Tp->savemode==1) break;
00336 }
00337 if (Tp->savemode==1) break;
00338 }
00339 }
00340
00341 delete [] r;
00342 }
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352 void linbartax::bf_matrix (matrix &n,double xi)
00353 {
00354 fillm (0.0,n);
00355 bf_lin_1d (n.a,xi);
00356 }
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 void linbartax::grad_matrix (matrix &gm,vector &x,double xi,double &jac)
00368 {
00369 long i;
00370 vector dx(nne);
00371 dx_bf_lin_1d (dx.a);
00372 derivatives_1d (dx,jac,x,xi);
00373
00374 for (i=0;i<nne;i++){
00375 gm[0][i]=dx[i];
00376 }
00377 }
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390 void linbartax::conductivity_matrix (long lcid,long eid,long ri,long ci,matrix &km)
00391 {
00392 long i,ii;
00393 double area,xinp,xi,ww,jac;
00394 ivector nodes(nne);
00395 vector x(nne),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),a(nne);
00396 matrix gm(ncomp,dofe[ri][ci]),d;
00397 matrix n(1,dofe[ri][ci]);
00398
00399
00400 Tt->give_elemnodes (eid,nodes);
00401
00402 Tc->give_area (eid,nodes,a);
00403
00404 Tt->give_node_coord1d (x,eid);
00405
00406 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00407 fillm (0.0,km);
00408
00409
00410 if (Tp->savemode==0)
00411 ii=Tt->elements[eid].ipp[ri][ci];
00412 if (Tp->savemode==1)
00413 ii=Tt->elements[eid].ipp[0][0];
00414
00415 for (i=0;i<intordkm[ri][ci];i++){
00416 xi=gp[i]; ww=w[i];
00417
00418
00419 grad_matrix (gm,x,xi,jac);
00420
00421
00422 reallocm(ncomp,ncomp,d);
00423 Tm->matcond (d,ii,ri,ci);
00424
00425 xinp = approx (xi,x);
00426
00427 area = approx (xi,a);
00428
00429 jac*=area*ww*xinp;
00430
00431
00432 bdbj (km.a,gm.a,d.a,jac,gm.m,gm.n);
00433
00434
00435 reallocm(1,ncomp,d);
00436
00437 Tm->matcond2(d,ii,ri,ci);
00438
00439 bf_matrix (n, xi);
00440 bdbjac(km, n, d, gm, jac);
00441
00442 ii++;
00443 }
00444
00445 if ((Tt->elements[eid].transi[lcid]==3) || (Tt->elements[eid].transi[lcid]==4)){
00446 transmission_matrix (lcid,eid,ri,ci,km);
00447 }
00448 }
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460 void linbartax::capacity_matrix (long eid,long ri,long ci,matrix &cm)
00461 {
00462 long i,ii;
00463 double jac,xi,ww,area,xinp,rho,c;
00464 ivector nodes(nne);
00465 vector x(nne),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),a(nne);
00466 matrix n(1,dofe[ri][ci]);
00467
00468
00469 Tt->give_elemnodes (eid,nodes);
00470
00471 Tc->give_area (eid,nodes,a);
00472
00473 Tt->give_node_coord1d (x,eid);
00474
00475 gauss_points (gp.a,w.a,intordkm[ri][ci]);
00476
00477 Tc->give_densitye (eid,rho);
00478
00479
00480 fillm (0.0,cm);
00481
00482 if (Tp->savemode==0)
00483 ii=Tt->elements[eid].ipp[ri][ci]+intordkm[ri][ci];
00484 if (Tp->savemode==1)
00485 ii=Tt->elements[eid].ipp[0][0]+intordkm[0][0];
00486
00487 for (i=0;i<intordcm[ri][ci];i++){
00488 xi=gp[i]; ww=w[i];
00489
00490 bf_matrix (n,xi);
00491
00492 c=Tm->capcoeff (ii,ri,ci);
00493 jac_1d (jac,x,xi);
00494 xinp = approx (xi,x);
00495
00496 area = approx (xi,a);
00497 jac*=area*ww*rho*c*xinp;;
00498
00499 nnj (cm.a,n.a,jac,n.m,n.n);
00500 ii++;
00501 }
00502
00503 }
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518 void linbartax::quantity_source_vector (vector &sv,vector &nodval,long eid,long ri,long ci)
00519 {
00520 long i;
00521 double jac,xi,ww,area,xinp;
00522 ivector nodes(nne);
00523 vector x(nne),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),v(dofe[ri][ci]),a(nne);
00524 matrix n(1,dofe[ri][ci]),nm(dofe[ri][ci],dofe[ri][ci]);
00525
00526
00527 Tt->give_elemnodes (eid,nodes);
00528
00529 Tc->give_area (eid,nodes,a);
00530
00531 Tt->give_node_coord1d (x,eid);
00532
00533 gauss_points (gp.a,w.a,intordcm[ri][ci]);
00534
00535 fillm (0.0,nm);
00536
00537 for (i=0;i<intordcm[ri][ci];i++){
00538 xi=gp[i]; ww=w[i];
00539
00540 bf_matrix (n,xi);
00541
00542 jac_1d (jac,x,xi);
00543
00544 xinp = approx (xi,x);
00545
00546
00547 area = approx (xi,a);
00548 jac*=ww*area*xinp;
00549
00550 nnj (nm.a,n.a,jac,n.m,n.n);
00551 }
00552 mxv (nm,nodval,v); addv (sv,v,sv);
00553
00554 }
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571 void linbartax::transmission_matrix (long lcid,long eid,long ri,long ci,matrix &km)
00572 {
00573 long i,leid,ii;
00574 double xi,area,xinp;
00575 double new_trc;
00576 long nn;
00577 bocontypet bc[2];
00578
00579 ivector nodes(nne);
00580 vector x(nne),trc(2),a(nne);
00581 matrix n(1,dofe[ri][ci]);
00582
00583
00584 Tt->give_elemnodes (eid,nodes);
00585
00586 Tc->give_area (eid,nodes,a);
00587
00588 Tt->give_node_coord1d (x,eid);
00589
00590 for (i=0;i<Tb->lc[lcid].neb;i++){
00591 if (Tb->lc[lcid].elemload[i].eid==eid){
00592 leid=i;
00593 break;
00594 }
00595 }
00596
00597
00598 Tb->lc[lcid].elemload[leid].give_bc (bc);
00599
00600 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
00601
00602
00603 if (Tp->savemode==0)
00604 ii=Tt->elements[eid].ipp[ri][ci];
00605 if (Tp->savemode==1)
00606 ii=Tt->elements[eid].ipp[0][0];
00607
00608
00609 if (bc[0]==4 || bc[0]>10){
00610 xi=-1.0;
00611
00612 bf_matrix (n,xi);
00613
00614
00615 nn = nodes[0];
00616
00617 Tm->transmission_transcoeff(new_trc,trc[0],ri,ci,nn,bc[0],ii);
00618
00619 trc[0]=new_trc;
00620
00621
00622 xinp = approx (xi,x);
00623
00624 area = approx (xi,a);
00625
00626 nnj (km.a,n.a,trc[0]*area*xinp,n.m,n.n);
00627 }
00628
00629 if (bc[1]==4 || bc[1]>10){
00630 xi=1.0;
00631
00632 bf_matrix (n,xi);
00633
00634
00635 nn = nodes[1];
00636
00637 Tm->transmission_transcoeff(new_trc,trc[1],ri,ci,nn,bc[1],ii);
00638
00639 trc[1]=new_trc;
00640
00641
00642 xinp = approx (xi,x);
00643
00644 area = approx (xi,a);
00645
00646 nnj (km.a,n.a,trc[1]*area*xinp,n.m,n.n);
00647 }
00648 }
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663 void linbartax::transmission_vector (vector &tmv,long lcid,long eid,long leid,long ri,long ci)
00664 {
00665 long ii;
00666 double xi,area,xinp;
00667 double new_nodval,tr;
00668 long nn;
00669 bocontypet bc[2];
00670
00671 ivector nodes(nne);
00672 vector x(nne),trc(2),trcn(nne),nodval(2),av(dofe[ri][ci]),v(dofe[ri][ci]),a(nne);
00673 vector trr(2);
00674 matrix n(1,dofe[ri][ci]),km(dofe[ri][ci],dofe[ri][ci]);
00675
00676
00677 Tt->give_elemnodes (eid,nodes);
00678
00679 Tc->give_area (eid,nodes,a);
00680
00681 Tt->give_node_coord1d (x,eid);
00682
00683
00684 Tb->lc[lcid].elemload[leid].give_bc (bc);
00685
00686 Tb->lc[lcid].elemload[leid].give_nodval (lcid,nodval);
00687
00688 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
00689
00690 Tb->lc[lcid].elemload[leid].give_trr (lcid,trr);
00691
00692
00693 if (Tp->savemode==0)
00694 ii=Tt->elements[eid].ipp[ri][ci];
00695 if (Tp->savemode==1)
00696 ii=Tt->elements[eid].ipp[0][0];
00697
00698 if (bc[0]==4 || bc[0]>10){
00699 fillm (0.0,km);
00700 xi=-1.0;
00701
00702 bf_matrix (n,xi);
00703
00704
00705 nn = nodes[0];
00706 tr = 0.0;
00707 if(bc[0] == 90)
00708 tr = trr[0]/trc[0];
00709
00710 Tm->transmission_nodval(new_nodval,nodval[0],tr,ri,ci,nn,bc[0],ii);
00711
00712 av[0]=new_nodval; av[1]=0.0;
00713
00714
00715 xinp = approx (xi,x);
00716
00717 area = approx (xi,a);
00718
00719 nnj (km.a,n.a,trc[0]*area*xinp,n.m,n.n);
00720
00721 mxv (km,av,v); addv (tmv,v,tmv);
00722 }
00723 if (bc[1]==4 || bc[1]>10){
00724 fillm (0.0,km);
00725 xi=1.0;
00726
00727 bf_matrix (n,xi);
00728
00729
00730 nn = nodes[1];
00731 tr = 0.0;
00732 if(bc[0] == 90)
00733 tr = trr[1]/trc[1];
00734
00735 Tm->transmission_nodval(new_nodval,nodval[1],tr,ri,ci,nn,bc[1],ii);
00736
00737 av[0]=0.0; av[1]=new_nodval;
00738
00739
00740 xinp = approx (xi,x);
00741
00742 area = approx (xi,a);
00743
00744 nnj (km.a,n.a,trc[1]*area*xinp,n.m,n.n);
00745
00746 mxv (km,av,v); addv (tmv,v,tmv);
00747 }
00748 }
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764 void linbartax::convection_vector (vector &f,long lcid,long eid,long leid,long ri,long ci)
00765 {
00766 bocontypet bc[2];
00767 double xi,area,xinp;
00768 ivector nodes(nne);
00769 vector x(nne),nodval(2),av(dofe[ri][ci]),v(dofe[ri][ci]),a(nne);
00770 matrix n(1,dofe[ri][ci]),km(dofe[ri][ci],dofe[ri][ci]);
00771
00772
00773 Tt->give_elemnodes (eid,nodes);
00774
00775 Tc->give_area (eid,nodes,a);
00776
00777 Tt->give_node_coord1d (x,eid);
00778
00779
00780 Tb->lc[lcid].elemload[leid].give_bc (bc);
00781
00782 Tb->lc[lcid].elemload[leid].give_nodval (lcid,nodval);
00783
00784
00785 if (bc[0]==2 || bc[0]==3 || bc[0]==4){
00786 fillm (0.0,km);
00787 xi=-1.0;
00788
00789 bf_matrix (n,xi);
00790
00791 xinp = approx (xi,x);
00792
00793 area = approx (xi,a);
00794
00795 nnj (km.a,n.a,area*xinp,n.m,n.n);
00796
00797 av[0]=nodval[0]; av[1]=0.0;
00798 mxv (km,av,v); addv (f,v,f);
00799 }
00800 if (bc[1]==2 || bc[1]==3 || bc[1]==4){
00801 fillm (0.0,km);
00802 xi=1.0;
00803 bf_matrix (n,xi);
00804
00805 xinp = approx (xi,x);
00806
00807 area = approx (xi,a);
00808
00809 nnj (km.a,n.a,area*xinp,n.m,n.n);
00810
00811 av[0]=0.0; av[1]=nodval[1];
00812 mxv (km,av,v); addv (f,v,f);
00813 }
00814 }
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 void linbartax::internal_fluxes (long lcid,long eid,vector &ifl)
00827 {
00828 long i,ipp;
00829 double jac,area,xinp,xi;
00830 ivector nodes(nne);
00831 vector x(nne),w(intordkm[lcid][lcid]),gp(intordkm[lcid][lcid]),fl(ncomp),contr(dofe[lcid][lcid]),a(nne);
00832 matrix gm(ncomp,dofe[lcid][lcid]),d(ncomp,ncomp);
00833
00834
00835 Tt->give_elemnodes (eid,nodes);
00836
00837 Tc->give_area (eid,nodes,a);
00838
00839 Tt->give_node_coord1d (x,eid);
00840
00841 gauss_points (gp.a,w.a,intordkm[lcid][lcid]);
00842
00843 if (Tp->savemode==0)
00844 ipp=Tt->elements[eid].ipp[lcid][lcid];
00845 if (Tp->savemode==1)
00846 ipp=Tt->elements[eid].ipp[0][0];
00847
00848
00849 for (i=0;i<intordkm[lcid][lcid];i++){
00850 xi=gp[i];
00851 xinp = approx (xi,x);
00852
00853 Tm->computenlfluxes (lcid,ipp);
00854
00855 Tm->givefluxes (lcid,ipp,fl);
00856
00857 grad_matrix (gm,x,gp[i],jac);
00858
00859 mtxv (gm,fl,contr);
00860
00861
00862 area = approx (xi,a);
00863 cmulv (area*jac*w[i]*xinp,contr);
00864
00865 addv (contr,ifl,ifl);
00866
00867 ipp++;
00868 }
00869 }
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880 void linbartax::res_conductivity_matrix (long eid,long lcid,matrix &km)
00881 {
00882 long i,j,*rcn,*ccn;
00883 matrix lkm;
00884
00885 for (i=0;i<ntm;i++){
00886 for (j=0;j<ntm;j++){
00887 rcn = new long [dofe[i][j]];
00888 ccn = new long [dofe[i][j]];
00889 reallocm (dofe[i][j],dofe[i][j],lkm);
00890 conductivity_matrix (i,eid,i,j,lkm);
00891 codnum (rcn,i);
00892 codnum (ccn,j);
00893 mat_localize (km,lkm,rcn,ccn);
00894 delete [] rcn; delete [] ccn;
00895 }
00896 }
00897 }
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907 void linbartax::res_capacity_matrix (long eid,matrix &cm)
00908 {
00909 long i,j,*rcn,*ccn,k,l;
00910 double s;
00911 matrix lcm;
00912
00913 for (i=0;i<ntm;i++){
00914 for (j=0;j<ntm;j++){
00915 rcn = new long [dofe[i][j]];
00916 ccn = new long [dofe[i][j]];
00917 reallocm (dofe[i][j],dofe[i][j],lcm);
00918 capacity_matrix (eid,i,j,lcm);
00919
00920
00921
00922 if (Tp->diagcap == 1){
00923 for (k=0;k<dofe[i][j];k++){
00924 s=0.0;
00925 for (l=0;l<dofe[i][j];l++){
00926 s+=lcm[k][l];
00927 lcm[k][l]=0.0;
00928 }
00929 lcm[k][k]=s;
00930 }
00931 }
00932
00933
00934 codnum (rcn,i);
00935 codnum (ccn,j);
00936 mat_localize (cm,lcm,rcn,ccn);
00937 delete [] rcn; delete [] ccn;
00938 }
00939 }
00940 }
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952 void linbartax::res_convection_vector (vector &f,long lcid,long eid,long leid)
00953 {
00954 long *cn;
00955 vector lf;
00956
00957 cn = new long [dofe[lcid][lcid]];
00958 reallocv (dofe[lcid][lcid],lf);
00959 convection_vector (lf,lcid,eid,leid,lcid,lcid);
00960 codnum (cn,lcid);
00961 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
00962 delete [] cn;
00963
00964 cmulv (-1.0,f,f);
00965
00966 }
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978 void linbartax::res_transmission_vector (vector &f,long lcid,long eid,long leid)
00979 {
00980 long *cn,i;
00981 vector lf;
00982
00983 cn = new long [dofe[lcid][lcid]];
00984 codnum (cn,lcid);
00985
00986
00987 for (i=0;i<ntm;i++){
00988 reallocv (dofe[lcid][lcid],lf);
00989
00990 if ((Tt->elements[eid].transi[i]==3) || (Tt->elements[eid].transi[i]==4))
00991 transmission_vector (lf,i,eid,leid,lcid,i);
00992
00993 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
00994 }
00995
00996 delete [] cn;
00997
00998 }
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009 void linbartax::res_quantity_source_vector (vector &sv,vector &nodval,long lcid,long eid)
01010 {
01011 long *cn;
01012 vector lsv;
01013
01014 cn = new long [dofe[lcid][lcid]];
01015 reallocv (dofe[lcid][lcid],lsv);
01016 codnum (cn,lcid);
01017
01018 quantity_source_vector (lsv,nodval,eid,lcid,lcid);
01019 locglob (sv.a,lsv.a,cn,dofe[lcid][lcid]);
01020
01021 delete [] cn;
01022 }
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032 void linbartax::res_internal_fluxes (long eid,vector &elemif)
01033 {
01034 long i,*cn;
01035 vector lif,tdnv;
01036 matrix cm;
01037
01038 for (i=0;i<ntm;i++){
01039 cn = new long [dofe[i][i]];
01040 reallocv (dofe[i][i],lif);
01041 codnum (cn,i);
01042 internal_fluxes (i,eid,lif);
01043 locglob (elemif.a,lif.a,cn,dofe[i][i]);
01044 delete [] cn;
01045 }
01046
01047 cn = new long [ndofe];
01048 Gtt->give_code_numbers (eid,cn);
01049 reallocv (ndofe,tdnv);
01050 reallocv (ndofe,lif);
01051 reallocm (ndofe,ndofe,cm);
01052
01053 nodalderivatives (tdnv.a,cn,ndofe);
01054 res_capacity_matrix (eid,cm);
01055 mxv (cm,tdnv,lif);
01056
01057 subv (elemif,lif,elemif);
01058
01059 }
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072 double linbartax::total_integral(long eid,vector &nodval)
01073 {
01074 long i;
01075 double area,xinp,xi,ww,jac,value,f;
01076 ivector nodes(nne);
01077 vector x(nne),w(2),gp(2),a(nne);
01078
01079
01080 Tt->give_elemnodes (eid,nodes);
01081
01082 Tc->give_area (eid,nodes,a);
01083
01084 Tt->give_node_coord1d (x,eid);
01085
01086 gauss_points (gp.a,w.a,2);
01087
01088 f = 0.0;
01089
01090 for (i=0;i<2;i++){
01091 xi=gp[i]; ww=w[i];
01092
01093 jac_1d (jac,x,xi);
01094
01095 xinp = approx (xi,x);
01096 value = approx (xi,nodval);
01097
01098 area = approx (xi,a);
01099
01100 jac*=area*ww*value*xinp;
01101
01102 f = f + jac;
01103 }
01104 return(f);
01105 }
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121 void linbartax::boundary_flux (vector &tmv,long lcid,long eid,long leid,long ri,long ci)
01122 {
01123 long ii;
01124 double xi,area,xinp;
01125 double new_nodval,tr;
01126 long nn;
01127 bocontypet bc[2];
01128
01129 ivector nodes(nne);
01130 vector x(nne),trc(2),trcn(nne),nodval(2),av(dofe[ri][ci]),v(dofe[ri][ci]),a(nne);
01131 vector trr(2);
01132 matrix n(1,dofe[ri][ci]),km(dofe[ri][ci],dofe[ri][ci]);
01133
01134
01135 Tt->give_elemnodes (eid,nodes);
01136
01137 Tc->give_area (eid,nodes,a);
01138
01139 Tt->give_node_coord1d (x,eid);
01140
01141
01142 Tb->lc[lcid].elemload[leid].give_bc (bc);
01143
01144 Tb->lc[lcid].elemload[leid].give_nodval (lcid,nodval);
01145
01146 Tb->lc[lcid].elemload[leid].give_trc (lcid,ci,trc);
01147
01148 Tb->lc[lcid].elemload[leid].give_trr (lcid,trr);
01149
01150 if (Tp->savemode==0)
01151 ii=Tt->elements[eid].ipp[ri][ci];
01152 if (Tp->savemode==1)
01153 ii=Tt->elements[eid].ipp[0][0];
01154
01155 if (bc[0]==4 || bc[0]>10){
01156 fillm (0.0,km);
01157 xi=-1.0;
01158
01159 bf_matrix (n,xi);
01160
01161 nn = nodes[0];
01162 tr = 0.0;
01163 if(bc[0] == 90)
01164 tr = trr[0]/trc[0];
01165
01166 Tm->transmission_flux(new_nodval,nodval[0],tr,ri,ci,nn,bc[0],ii);
01167
01168 av[0]=new_nodval; av[1]=0.0;
01169
01170 xinp = approx (xi,x);
01171
01172 area = approx (xi,a);
01173
01174 nnj (km.a,n.a,trc[0]*area*xinp,n.m,n.n);
01175
01176 mxv (km,av,v); addv (tmv,v,tmv);
01177 }
01178 if (bc[1]==4 || bc[1]>10){
01179 fillm (0.0,km);
01180 xi=1.0;
01181
01182 bf_matrix (n,xi);
01183
01184 nn = nodes[1];
01185 tr = 0.0;
01186 if(bc[0] == 90)
01187 tr = trr[1]/trc[1];
01188
01189 Tm->transmission_flux(new_nodval,nodval[1],tr,ri,ci,nn,bc[1],ii);
01190
01191 av[0]=0.0; av[1]=new_nodval;
01192
01193 xinp = approx (xi,x);
01194
01195 area = approx (xi,a);
01196
01197 nnj (km.a,n.a,trc[1]*area*xinp,n.m,n.n);
01198
01199 mxv (km,av,v); addv (tmv,v,tmv);
01200 }
01201 }
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215 void linbartax::res_boundary_flux (vector &f,long lcid,long eid,long leid)
01216 {
01217 long *cn,i;
01218 vector lf;
01219
01220 cn = new long [dofe[lcid][lcid]];
01221 codnum (cn,lcid);
01222
01223
01224 for (i=0;i<ntm;i++){
01225 reallocv (dofe[lcid][lcid],lf);
01226
01227 if ((Tt->elements[eid].transi[i]==3) || (Tt->elements[eid].transi[i]==4))
01228 boundary_flux (lf,i,eid,leid,lcid,i);
01229
01230 locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
01231 }
01232
01233 delete [] cn;
01234 }
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 void linbartax::nod_others (long lcid,long eid,long ri,long ci)
01247 {
01248 long i,j,k,ipp,ncompother,ii;
01249 vector other,h,r(ndofe);
01250 ivector nod(nne),cn(ndofe);
01251
01252 Tt->give_elemnodes (eid,nod);
01253 Tt->give_code_numbers (eid,cn.a);
01254 elemvalues (lcid,eid,r.a,cn.a,ndofe);
01255 ncompother = Tm->givencompother();
01256 reallocv (ncompother,other);
01257 reallocv (ntm,h);
01258
01259 ii = 0;
01260 for (i=0;i<nne;i++){
01261 if (Tp->savemode==0)
01262 ipp=Tt->elements[eid].ipp[ri][ci];
01263 if (Tp->savemode==1)
01264 ipp=Tt->elements[eid].ipp[0][0];
01265
01266 for (j=0;j<ntm;j++){
01267 h[j] = r[ii];
01268 ii++;
01269 }
01270
01271 for(k=0;k<ncompother;k++)
01272 other[k] = Tm->givecompother(k,ipp,h.a);
01273
01274 }
01275 }