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