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