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