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