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