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