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