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