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