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