00001 #include "transtop.h"
00002 #include "globalt.h"
00003 #include "globmatt.h"
00004 #include "elemswitcht.h"
00005 #include <string.h>
00006 #include <limits.h>
00007 #include <stdio.h>
00008
00009 transtop::transtop ()
00010 {
00011
00012 nn=0;
00013
00014 ne=0;
00015
00016
00017 nodes=NULL;
00018
00019 elements=NULL;
00020
00021 enodes=NULL;
00022
00023 edges=NULL;
00024
00025
00026
00027 nedis=0;
00028 lnnd=NULL;
00029 lnd=NULL;
00030
00031
00032 ncycl = NULL;
00033
00034
00035 mnn = 0;
00036
00037
00038
00039
00040
00041 }
00042
00043 transtop::~transtop ()
00044 {
00045 long i;
00046
00047 if (nodes!=NULL)
00048 delete [] nodes;
00049 if (elements!=NULL)
00050 delete [] elements;
00051 if (enodes!=NULL)
00052 delete [] enodes;
00053 if (edges!=NULL)
00054 delete [] edges;
00055 if (lnnd!=NULL)
00056 delete [] lnnd;
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 if (lnd!=NULL){
00071 for (i=0;i<nedis;i++){
00072 delete [] lnd[i];
00073 }
00074 delete [] lnd;
00075 }
00076 if (ncycl!=NULL)
00077 delete [] ncycl;
00078 }
00079
00080
00081
00082
00083
00084
00085 void transtop::read (XFILE *in)
00086 {
00087 long i,j,k,ndofn;
00088
00089
00090
00091
00092 xfscanf (in,"%k%ld","number_of_nodes",&nn);
00093 Gtt->alloc_nodes (nn);
00094
00095 if (Mesprt==1) fprintf (stdout,"\n number of nodes %ld",nn);
00096 nodes = new nodet [nn];
00097 for (i=0;i<nn;i++){
00098 j=nn+1;
00099 xfscanf (in,"%ld",&j);
00100 if (j<1)
00101 print_err("node number in function read is less than 1",__FILE__,__LINE__,__func__);
00102 if (j>nn)
00103 print_err("node number in function read is greater than number of nodes",__FILE__,__LINE__,__func__);
00104 Gtt->gnodes[j-1].read (in);
00105 ndofn = Gtt->gnodes[j-1].ndofn;
00106 nodes[j-1].read (in,ndofn);
00107 }
00108
00109
00110
00111
00112 xfscanf (in,"%k%ld","number_of_constraints",&ncn);
00113 if (Mesprt==1) fprintf (stdout,"\n number of constrained nodes %ld",ncn);
00114 for (i=0;i<ncn;i++){
00115 k=nn+1;
00116 xfscanf (in,"%ld",&k);
00117 if (k<1)
00118 print_err("number of constrained node in function read is less than 1",__FILE__,__LINE__,__func__);
00119 if (k>nn){
00120 print_err("number of constrained node in function read is greater than number of all nodes",__FILE__,__LINE__,__func__);
00121 }
00122 Gtt->gnodes[k-1].constr (Gtt->dofcontr,in);
00123 }
00124
00125
00126
00127
00128
00129 if(Tp->homogt == 1){
00130 xfscanf (in,"%k%ld","masternode_number",&mnn);
00131 if (Mesprt==1) fprintf (stdout,"\n master node number for homogenization %ld",mnn);
00132 mnn = mnn-1;
00133 }
00134
00135
00136
00137
00138 xfscanf (in,"%k%ld","number_of_elements",&ne);
00139 Gtt->alloc_elements (ne);
00140 if (Mesprt==1) fprintf (stdout,"\n number of elements %ld",ne);
00141 elements = new elementt [ne];
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 for (i=0;i<ne;i++){
00164 j=ne+1;
00165 xfscanf (in,"%ld",&j);
00166 if (j<1)
00167 print_err("element number in function read is less than 1",__FILE__,__LINE__,__func__);
00168 if (j>ne){
00169 print_err("element number in function read is greater than total number of elements",__FILE__,__LINE__,__func__);
00170 }
00171 elements[j-1].read (in,j-1);
00172 }
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 Gtt->lneso_init ();
00189
00190
00191 if (Tp->tprob == growing_np_problem){
00192 Gtt->auxinf_init ();
00193 }
00194 if (Tp->tprob == growing_np_problem_nonlin){
00195 Gtt->auxinf_init ();
00196 }
00197
00198
00199 if (Gtt->rst==1){
00200
00201
00202
00203 Gtt->read_seq_top (in);
00204 }
00205 if (Gtt->rst==2){
00206
00207
00208 Tp->ssle->feti.read_booldata (in);
00209 }
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265 }
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 void transtop::print (FILE *out)
00276 {
00277 long i;
00278
00279
00280
00281
00282 fprintf (out,"\n%ld #number of nodes\n\n",nn);
00283
00284 for (i=0;i<nn;i++){
00285 fprintf (out,"\n %ld",i+1);
00286 Gtt->gnodes[i].print (out);
00287 nodes[i].print (out);
00288 }
00289
00290
00291
00292
00293 fprintf (out,"\n%ld #number of constrain nodes\n\n",ncn);
00294 for (i=0;i<ncn;i++){
00295 fprintf (out,"\n %ld",i);
00296 Gtt->gnodes[i].print_constr (Gtt->dofcontr,out);
00297 }
00298
00299
00300
00301
00302
00303 if(Tp->homogt == 1){
00304 fprintf (out,"\n %ld #master node number\n\n",mnn);
00305 }
00306
00307
00308
00309
00310 fprintf (out,"\n %ld #number of elements\n\n",ne);
00311 for (i=0;i<ne;i++){
00312 fprintf (out,"\n %ld",i+1);
00313
00314 elements[i].print (out,i);
00315 }
00316 }
00317
00318
00319
00320
00321
00322
00323
00324 elemtypet transtop::give_elem_type (long eid)
00325 {
00326 return (elements[eid].te);
00327 }
00328
00329
00330
00331
00332
00333
00334
00335 long transtop::give_dof (long nid,long n)
00336 {
00337 return Gtt->give_dof (nid,n);
00338 }
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 void transtop::give_elemnodes (long eid,ivector &nodes)
00349 {
00350 Gtt->give_nodes (eid,nodes);
00351 }
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361 void transtop::give_code_numbers (long eid,long *cn)
00362 {
00363 Gtt->give_code_numbers (eid,cn);
00364 }
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375 void transtop::give_medium_code_numbers (long eid,long ri,long *cn)
00376 {
00377 long i,ndofe,nne,ntm;
00378 long *ecn;
00379
00380
00381 ndofe = give_ndofe (eid);
00382 ecn = new long [ndofe];
00383
00384 Gtt->give_code_numbers (eid,ecn);
00385
00386
00387 nne = give_nne (eid);
00388
00389
00390 ntm = Tp->ntm;
00391
00392 for (i=0;i<nne;i++){
00393 cn[i]=ecn[i*ntm+ri];
00394 }
00395
00396 delete [] ecn;
00397 }
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 void transtop::give_node_code_numbers (long nid,long *cn)
00408 {
00409 Gtt->give_node_code_numbers (nid,cn);
00410 }
00411
00412
00413
00414
00415
00416
00417 void transtop::give_node_coord1d (vector &x,long eid)
00418 {
00419 Gtt->give_node_coord1d (x,eid);
00420 }
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430 void transtop::give_node_coord2d (vector &x,vector &y,long eid)
00431 {
00432 Gtt->give_node_coord2d (x,y,eid);
00433 }
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443 void transtop::give_node_coord3d (vector &x,vector &y,vector &z,long eid)
00444 {
00445 Gtt->give_node_coord3d (x,y,z,eid);
00446 }
00447
00448
00449
00450
00451
00452
00453
00454
00455 long transtop::give_ndofe (long eid)
00456 {
00457 long ndofe;
00458 elemtypet te;
00459
00460 te = give_elem_type (eid);
00461
00462 switch (te){
00463 case barlint:{ ndofe=Lbt->ndofe; break; }
00464 case barlintax:{ ndofe=Lbat->ndofe; break; }
00465 case barquadt:{ ndofe=Qbt->ndofe; break; }
00466 case barquadtax:{ ndofe=Qbat->ndofe; break; }
00467 case trlint:{ ndofe=Ltt->ndofe; break; }
00468 case trlaxisym:{ ndofe=Ltat->ndofe; break; }
00469 case quadlint:{ ndofe=Lqt->ndofe; break; }
00470 case quadquadt:{ ndofe=Qqt->ndofe; break; }
00471 case quadquadtax:{ ndofe=Qqat->ndofe; break; }
00472 case quadlaxisym:{ ndofe=Lqat->ndofe; break; }
00473 case lineartett:{ ndofe=Ltett->ndofe; break; }
00474 case linearhext:{ ndofe=Lht->ndofe; break; }
00475 case quadratichext:{ ndofe=Qht->ndofe; break; }
00476 case gen2del:{ ndofe=G2d->ndofe; break; }
00477 default:{
00478 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00479 }
00480 }
00481 return ndofe;
00482
00483 }
00484
00485
00486
00487
00488
00489
00490
00491 long transtop::give_ndofn (long nid)
00492 {
00493 return Gtt->give_ndofn (nid);
00494 }
00495
00496
00497
00498
00499
00500
00501
00502
00503 long transtop::give_tnip (long eid)
00504 {
00505 if (eid < ne-1)
00506 return elements[eid+1].ipp[0][0] - elements[eid].ipp[0][0];
00507
00508 return Tm->tnip-elements[eid].ipp[0][0];
00509 }
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520 long transtop::give_ncomp (long eid)
00521 {
00522 long ncomp;
00523 elemtypet te;
00524
00525 te = give_elem_type (eid);
00526
00527 switch (te){
00528 case barlint:{ ncomp=Lbt->ncomp; break; }
00529 case barlintax:{ ncomp=Lbat->ncomp; break; }
00530 case barquadt:{ ncomp=Qbt->ncomp; break; }
00531 case barquadtax:{ ncomp=Qbat->ncomp; break; }
00532 case trlint:{ ncomp=Ltt->ncomp; break; }
00533 case trlaxisym:{ ncomp=Ltat->ncomp; break; }
00534 case quadlint:{ ncomp=Lqt->ncomp; break; }
00535 case quadquadt:{ ncomp=Qqt->ncomp; break; }
00536 case quadquadtax:{ ncomp=Qqat->ncomp; break; }
00537 case quadlaxisym:{ ncomp=Lqat->ncomp; break; }
00538 case lineartett:{ ncomp=Ltett->ncomp; break; }
00539 case linearhext:{ ncomp=Lht->ncomp; break; }
00540 case quadratichext:{ ncomp=Qht->ncomp; break; }
00541 case gen2del:{ ncomp=G2d->ncomp; break; }
00542
00543 default:{
00544 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00545 }
00546 }
00547 return ncomp;
00548 }
00549
00550
00551
00552
00553
00554
00555
00556
00557 long transtop::give_ned (long eid)
00558 {
00559 long ned;
00560 elemtypet te;
00561
00562 te = give_elem_type (eid);
00563
00564 switch (te){
00565 case barlint:{ ned=Lbt->ned; break; }
00566 case barlintax:{ ned=Lbat->ned; break; }
00567 case barquadt:{ ned=Qbt->ned; break; }
00568 case barquadtax:{ ned=Qbat->ned; break; }
00569 case trlint:{ ned=Ltt->ned; break; }
00570 case trlaxisym:{ ned=Ltat->ned; break; }
00571 case quadlint:{ ned=Lqt->ned; break; }
00572 case quadlaxisym:{ ned=Lqat->ned; break; }
00573 case quadquadt:{ ned=Qqt->ned; break; }
00574 case quadquadtax:{ ned=Qqat->ned; break; }
00575 case lineartett:{ ned=Ltett->ned; break; }
00576 case linearhext:{ ned=Lht->ned; break; }
00577 case quadratichext:{ned=Qht->ned; break; }
00578 case gen2del:{ ned=0; break; }
00579
00580 default:{
00581 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00582 }
00583 }
00584 return ned;
00585
00586 }
00587
00588
00589
00590
00591
00592
00593
00594
00595 long transtop::give_nned (long eid)
00596 {
00597 long nned;
00598 elemtypet te;
00599
00600 te = give_elem_type (eid);
00601
00602 switch (te){
00603 case barlint:{ nned=Lbt->nned; break; }
00604 case barlintax:{ nned=Lbat->nned; break; }
00605 case barquadt:{ nned=Qbt->nned; break; }
00606 case barquadtax:{ nned=Qbat->nned; break; }
00607 case trlint:{ nned=Ltt->nned; break; }
00608 case trlaxisym:{ nned=Ltat->nned; break; }
00609 case quadlint:{ nned=Lqt->nned; break; }
00610 case quadquadt:{ nned=Qqt->nned; break; }
00611 case quadquadtax:{ nned=Qqat->nned; break; }
00612 case quadlaxisym:{ nned=Lqat->nned; break; }
00613 case lineartett:{ nned=Ltett->nned; break; }
00614 case linearhext:{ nned=Lht->nned; break; }
00615 case quadratichext:{ nned=Qht->nned; break; }
00616 case gen2del:{ nned=0; break; }
00617 default:{
00618 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00619 }
00620 }
00621 return nned;
00622 }
00623
00624
00625
00626
00627
00628
00629
00630
00631 long transtop::give_nsurf (long eid)
00632 {
00633 long nsurf;
00634 elemtypet te;
00635
00636 te = give_elem_type (eid);
00637
00638 switch (te){
00639 case barlint:{ nsurf=0; break; }
00640 case barlintax:{ nsurf=0; break; }
00641 case barquadt:{ nsurf=0; break; }
00642 case barquadtax:{ nsurf=0; break; }
00643 case trlint:{ nsurf=0; break; }
00644 case trlaxisym:{ nsurf=0; break; }
00645 case quadlint:{ nsurf=0; break; }
00646 case quadquadt:{ nsurf=0; break; }
00647 case quadquadtax:{ nsurf=0; break; }
00648 case quadlaxisym:{ nsurf=0; break; }
00649 case gen2del:{ nsurf=0; break; }
00650
00651 case lineartett:{ nsurf=Ltett->nsurf; break; }
00652 case linearhext:{ nsurf=Lht->nsurf; break; }
00653 case quadratichext:{ nsurf=Qht->nsurf; break; }
00654 default:{
00655 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00656 }
00657 }
00658 return nsurf;
00659
00660 }
00661
00662
00663
00664
00665
00666
00667
00668
00669 long transtop::give_nnsurf (long eid)
00670 {
00671 long nnsurf;
00672 elemtypet te;
00673
00674 te = give_elem_type (eid);
00675
00676 switch (te){
00677 case lineartett:{ nnsurf=Ltett->nnsurf; break; }
00678 case linearhext:{ nnsurf=Lht->nnsurf; break; }
00679 case quadratichext:{ nnsurf=Qht->nnsurf; break; }
00680 default:{
00681 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00682 }
00683 }
00684 return nnsurf;
00685
00686 }
00687
00688
00689
00690
00691
00692
00693
00694
00695 long transtop::give_nne (long eid)
00696 {
00697 long nne;
00698 elemtypet te;
00699
00700 te = give_elem_type (eid);
00701
00702 switch (te){
00703 case barlint:{ nne=Lbt->nne; break; }
00704 case barlintax:{ nne=Lbat->nne; break; }
00705 case barquadt:{ nne=Qbt->nne; break; }
00706 case barquadtax:{ nne=Qbat->nne; break; }
00707 case trlint:{ nne=Ltt->nne; break; }
00708 case trlaxisym:{ nne=Ltat->nne; break; }
00709 case quadlint:{ nne=Lqt->nne; break; }
00710 case quadquadt:{ nne=Qqt->nne; break; }
00711 case quadquadtax:{ nne=Qqat->nne; break; }
00712 case quadlaxisym:{ nne=Lqat->nne; break; }
00713 case lineartett:{ nne=Ltett->nne; break; }
00714 case linearhext:{ nne=Lht->nne; break; }
00715 case quadratichext:{ nne=Qht->nne; break; }
00716 case gen2del:{ nne=G2d->nne; break; }
00717 default:{
00718 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00719 }
00720 }
00721 return nne;
00722
00723 }
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734 long transtop::give_nip (long eid,long ri,long ci)
00735 {
00736 long nip;
00737 elemtypet te;
00738
00739 te = give_elem_type (eid);
00740
00741 switch (te){
00742 case barlint:{ nip=Lbt->nip[ri][ci]; break; }
00743 case barlintax:{ nip=Lbat->nip[ri][ci]; break; }
00744 case barquadt:{ nip=Qbt->nip[ri][ci]; break; }
00745 case barquadtax:{ nip=Qbat->nip[ri][ci]; break; }
00746 case trlint:{ nip=Ltt->nip[ri][ci]; break; }
00747 case trlaxisym:{ nip=Ltat->nip[ri][ci]; break; }
00748 case quadlint:{ nip=Lqt->nip[ri][ci]; break; }
00749 case quadquadt:{ nip=Qqt->nip[ri][ci]; break; }
00750 case quadquadtax:{ nip=Qqat->nip[ri][ci]; break; }
00751 case quadlaxisym:{ nip=Lqat->nip[ri][ci]; break; }
00752 case lineartett:{ nip=Ltett->nip[ri][ci]; break; }
00753 case linearhext:{ nip=Lht->nip[ri][ci]; break; }
00754 case quadratichext:{ nip=Qht->nip[ri][ci]; break; }
00755 case gen2del:{ nip=G2d->nip[ri][ci]; break; }
00756 default:{
00757 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00758 }
00759 }
00760 return nip;
00761
00762 }
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776 long transtop::give_intordkm(long eid, long ri, long ci)
00777 {
00778 elemtypet te;
00779 long ret = 0;
00780
00781 te=Tt->give_elem_type (eid);
00782
00783 switch (te){
00784 case barlint:{
00785 ret = Lbt->intordkm[ri][ci];
00786 break;
00787 }
00788 case barlintax:{
00789 ret = Lbat->intordkm[ri][ci];
00790 break;
00791 }
00792 case barquadt:{
00793 ret = Qbt->intordkm[ri][ci];
00794 break;
00795 }
00796 case barquadtax:{
00797 ret = Qbat->intordkm[ri][ci];
00798 break;
00799 }
00800 case trlint:{
00801 ret = Ltt->intordkm[ri][ci];
00802 break;
00803 }
00804 case trlaxisym:{
00805 ret = Ltat->intordkm[ri][ci];
00806 break;
00807 }
00808 case quadlint:{
00809 ret = Lqt->intordkm[ri][ci];
00810 break;
00811 }
00812 case quadquadt:{
00813 ret = Qqt->intordkm[ri][ci];
00814 break;
00815 }
00816 case quadquadtax:{
00817 ret = Qqat->intordkm[ri][ci];
00818 break;
00819 }
00820 case quadlaxisym:{
00821 ret = Lqat->intordkm[ri][ci];
00822 break;
00823 }
00824 case lineartett:{
00825 ret = Ltett->intordkm[ri][ci];
00826 break;
00827 }
00828 case linearhext:{
00829 ret = Lht->intordkm[ri][ci];
00830 break;
00831 }
00832 case quadratichext:{
00833 ret = Qht->intordkm[ri][ci];
00834 break;
00835 }
00836 default:{
00837 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00838 }
00839 }
00840 return ret;
00841 }
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857 long transtop::give_intordcm(long eid, long ri, long ci)
00858 {
00859 elemtypet te;
00860 long ret = 0;
00861
00862 te=Tt->give_elem_type (eid);
00863
00864 switch (te){
00865 case barlint:{
00866 ret = Lbt->intordcm[ri][ci];
00867 break;
00868 }
00869 case barlintax:{
00870 ret = Lbat->intordcm[ri][ci];
00871 break;
00872 }
00873 case barquadt:{
00874 ret = Qbt->intordcm[ri][ci];
00875 break;
00876 }
00877 case barquadtax:{
00878 ret = Qbat->intordcm[ri][ci];
00879 break;
00880 }
00881 case trlint:{
00882 ret = Ltt->intordcm[ri][ci];
00883 break;
00884 }
00885 case trlaxisym:{
00886 ret = Ltat->intordcm[ri][ci];
00887 break;
00888 }
00889 case quadlint:{
00890 ret = Lqt->intordcm[ri][ci];
00891 break;
00892 }
00893 case quadquadt:{
00894 ret = Qqt->intordcm[ri][ci];
00895 break;
00896 }
00897 case quadquadtax:{
00898 ret = Qqat->intordcm[ri][ci];
00899 break;
00900 }
00901 case quadlaxisym:{
00902 ret = Lqat->intordcm[ri][ci];
00903 break;
00904 }
00905 case lineartett:{
00906 ret = Ltett->intordcm[ri][ci];
00907 break;
00908 }
00909 case linearhext:{
00910 ret = Lht->intordcm[ri][ci];
00911 break;
00912 }
00913 case quadratichext:{
00914 ret = Qht->intordcm[ri][ci];
00915 break;
00916 }
00917 default:{
00918 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00919 }
00920 }
00921 return ret;
00922 }
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934 long transtop::give_degree (long eid)
00935 {
00936 long deg;
00937 elemtypet te;
00938
00939 te = give_elem_type (eid);
00940
00941 switch (te){
00942 case barlint:{ deg=1; break; }
00943 case barlintax:{ deg=1; break; }
00944 case barquadt:{ deg=2; break; }
00945 case barquadtax:{ deg=2; break; }
00946 case trlint:{ deg=1; break; }
00947 case trlaxisym:{ deg=1; break; }
00948 case quadlint:{ deg=1; break; }
00949 case quadquadt:{ deg=2; break; }
00950 case quadquadtax:{ deg=2; break; }
00951 case quadlaxisym:{ deg=1; break; }
00952 case lineartett:{ deg=1; break; }
00953 case linearhext:{ deg=1; break; }
00954 case quadratichext:{ deg=2; break; }
00955 default:{
00956 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00957 }
00958 }
00959 return deg;
00960 }
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970 long transtop::give_dimension (long eid)
00971 {
00972 long dim;
00973 elemtypet te;
00974
00975 te = give_elem_type (eid);
00976
00977 switch (te){
00978 case barlint:{ dim=1; break; }
00979 case barlintax:{ dim=1; break; }
00980 case barquadt:{ dim=1; break; }
00981 case barquadtax:{ dim=1; break; }
00982 case trlint:{ dim=2; break; }
00983 case trlaxisym:{ dim=2; break; }
00984 case quadlint:{ dim=2; break; }
00985 case quadquadt:{ dim=2; break; }
00986 case quadquadtax:{ dim=2; break; }
00987 case quadlaxisym:{ dim=2; break; }
00988 case lineartett:{ dim=3; break; }
00989 case linearhext:{ dim=3; break; }
00990 case quadratichext:{ dim=3; break; }
00991 case gen2del:{ dim=2; break; }
00992 default:{
00993 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
00994 }
00995 }
00996 return dim;
00997 }
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009 void transtop::give_end_nodes (long eid,long enid,long *nodes)
01010 {
01011 long nne;
01012 elemtypet te;
01013 ivector nod,endnod;
01014
01015
01016 te = give_elem_type (eid);
01017
01018 nne = give_nne (eid);
01019
01020 allocv (nne,nod);
01021 allocv (1,endnod);
01022
01023
01024 give_elemnodes (eid,nod);
01025
01026 switch (te){
01027 case barlint:{
01028 endnod[0]=enid;
01029 break;
01030 }
01031 case barlintax:{
01032 endnod[0]=enid;
01033 break;
01034 }
01035 case barquadt:{
01036 endnod[0]=enid;
01037 break;
01038 }
01039 case barquadtax:{
01040 endnod[0]=enid;
01041 break;
01042 }
01043 case trlint:{
01044 break;
01045 }
01046 case trlaxisym:{
01047 break;
01048 }
01049 case quadlint:{
01050 break;
01051 }
01052 case quadlaxisym:{
01053 break;
01054 }
01055 case quadquadt:{
01056 break;
01057 }
01058 case lineartett:{
01059 break;
01060 }
01061 case linearhext:{
01062 break;
01063 }
01064 case quadratichext:{
01065 break;
01066 }
01067 case gen2del:{
01068 break;
01069 }
01070 default:{
01071 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01072 }
01073 }
01074
01075 nodes[0]=nod[endnod[0]];
01076
01077 destrv (nod); destrv (endnod);
01078 }
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090 void transtop::give_edge_nodes (long eid,long edid,long *nodes)
01091 {
01092 long i,nne,nned;
01093 elemtypet te;
01094 ivector nod,edgenod;
01095
01096
01097 te = give_elem_type (eid);
01098
01099 nne = give_nne (eid);
01100
01101 nned = give_nned (eid);
01102
01103 allocv (nne,nod);
01104 allocv (nned,edgenod);
01105
01106
01107 give_elemnodes (eid,nod);
01108
01109 switch (te){
01110 case barlint:{
01111 edgenod[0]=edid;
01112 break;
01113 }
01114 case barlintax:{
01115 break;
01116 }
01117 case barquadt:{
01118 break;
01119 }
01120 case barquadtax:{
01121 break;
01122 }
01123 case trlint:{
01124 lintriangle_edgnod (edgenod.a,edid);
01125 break;
01126 }
01127 case trlaxisym:{
01128 lintriangle_edgnod (edgenod.a,edid);
01129 break;
01130 }
01131 case quadlint:{
01132 linquadrilat_edgnod (edgenod.a,edid);
01133 break;
01134 }
01135 case quadlaxisym:{
01136 linquadrilat_edgnod (edgenod.a,edid);
01137 break;
01138 }
01139 case quadquadt:{
01140 quadquadrilat_edgnod (edgenod.a,edid);
01141 break;
01142 }
01143 case lineartett:{
01144 lintetrahedral_edgnod(edgenod.a,edid);
01145 break;
01146 }
01147 case linearhext:{
01148 linhexahedral_edgnod(edgenod.a,edid);
01149 break;
01150 }
01151 case quadratichext:{
01152 quadhexahedral_edgnod (edgenod.a,edid);
01153 break;
01154 }
01155 case gen2del:{
01156 break;
01157 }
01158 default:{
01159 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01160 }
01161 }
01162
01163 for (i=0;i<nned;i++){
01164 nodes[i]=nod[edgenod[i]];
01165 }
01166 }
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177 void transtop::give_surface_nodes (long eid,long surfid,long *nodes)
01178 {
01179 long i,nne,nnsurf;
01180 elemtypet te;
01181 ivector nod,surfnod;
01182
01183
01184 te = give_elem_type (eid);
01185
01186 nne = give_nne (eid);
01187
01188 nnsurf = give_nnsurf (eid);
01189
01190 allocv (nne,nod);
01191 allocv (nnsurf,surfnod);
01192
01193
01194 give_elemnodes (eid,nod);
01195
01196 switch (te){
01197 case barlint:{
01198 break;
01199 }
01200 case barlintax:{
01201 break;
01202 }
01203 case barquadt:{
01204 break;
01205 }
01206 case barquadtax:{
01207 break;
01208 }
01209 case trlint:{
01210 break;
01211 }
01212 case trlaxisym:{
01213 break;
01214 }
01215 case quadlint:{
01216 break;
01217 }
01218 case quadlaxisym:{
01219 break;
01220 }
01221
01222 case lineartett:{
01223 break;
01224 }
01225 case linearhext:{
01226 linhexahedral_surfnod (surfnod.a,surfid);
01227 break;
01228 }
01229 case quadratichext:{
01230 quadhexahedral_surfnod (surfnod.a,surfid);
01231 break;
01232 }
01233 case gen2del:{
01234 break;
01235 }
01236 default:{
01237 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01238 }
01239 }
01240
01241 for (i=0;i<nnsurf;i++){
01242 nodes[i]=nod[surfnod[i]];
01243 }
01244
01245 destrv (nod); destrv (surfnod);
01246 }
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258 void transtop::give_nbobjects (long eid,long &bco,long &ncompbco)
01259 {
01260
01261 switch (elements[eid].te){
01262
01263 case barlint:{
01264 bco=Lbt->nen;
01265 ncompbco=1;
01266 break;
01267 }
01268 case barlintax:{
01269 bco=Lbat->nen;
01270 ncompbco=1;
01271 break;
01272 }
01273 case barquadt:{
01274 bco=Qbt->nen;
01275 ncompbco=1;
01276 break;
01277 }
01278 case barquadtax:{
01279 bco=Qbat->nen;
01280 ncompbco=1;
01281 break;
01282 }
01283
01284
01285 case trlint:{
01286 bco=Ltt->ned;
01287 ncompbco=Ltt->nned;
01288 break;
01289 }
01290 case trlaxisym:{
01291 bco=Ltat->ned;
01292 ncompbco=Ltat->nned;
01293 break;
01294 }
01295 case quadlint:{
01296 bco=Lqt->ned;
01297 ncompbco=Lqt->nned;
01298 break;
01299 }
01300 case quadlaxisym:{
01301 bco=Lqat->ned;
01302 ncompbco=Lqat->nned;
01303 break;
01304 }
01305 case quadquadt:{
01306 bco=Qqt->ned;
01307 ncompbco=Qqt->nned;
01308 break;
01309 }
01310 case quadquadtax:{
01311 bco=Qqat->ned;
01312 ncompbco=Qqat->nned;
01313 break;
01314 }
01315
01316
01317 case lineartett:{
01318 bco=Ltett->nsurf;
01319 ncompbco=Ltett->nnsurf;
01320 break;
01321 }
01322 case linearhext:{
01323 bco=Lht->nsurf;
01324 ncompbco=Lht->nnsurf;
01325 break;
01326 }
01327 case quadratichext:{
01328 bco=Qht->nsurf;
01329 ncompbco=Qht->nnsurf;
01330 break;
01331 }
01332 case gen2del:{
01333 bco=0;
01334 ncompbco=0;
01335 break;
01336 }
01337 default:{
01338 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01339 }
01340 }
01341
01342 }
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353 void transtop::give_bonodes (long eid,long boid,long *nod)
01354 {
01355
01356 switch (elements[eid].te){
01357
01358 case barlint:
01359 case barlintax:
01360 case barquadt:
01361 case barquadtax:{
01362 give_end_nodes (eid,boid,nod);
01363 break;
01364 }
01365
01366
01367 case trlint:
01368 case trlaxisym:
01369 case quadlint:
01370 case quadlaxisym:
01371 case quadquadt:
01372 case quadquadtax:{
01373 give_edge_nodes (eid,boid,nod);
01374 break;
01375 }
01376
01377
01378 case lineartett:
01379 case linearhext:
01380 case quadratichext:{
01381 give_surface_nodes (eid,boid,nod);
01382 break;
01383 }
01384 case gen2del:{
01385 break;
01386 }
01387 default:{
01388 print_err("unknown element type is required",__FILE__,__LINE__,__func__);
01389 }
01390 }
01391
01392 }
01393
01394
01395
01396
01397 void transtop::alloc_prep(long nn, long ne)
01398 {
01399 nodes = new nodet[nn];
01400 elements = new elementt[ne];
01401 }
01402
01403
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
01436
01437
01438 void transtop::alloc_nodes (void)
01439 {
01440 long i,j,nne,ncomp, ncompo, ncompeqo, ipp;
01441 ivector nod;
01442
01443
01444 for (i=0;i<ne;i++){
01445
01446 nne = give_nne (i);
01447 allocv (nne,nod);
01448
01449 give_elemnodes (i,nod);
01450
01451 ncomp=give_dimension (i);
01452
01453 ncompo=Tm->givencompother ();
01454
01455
01456 ipp = elements[i].ipp[0][0];
01457
01458 ncompeqo=Tm->givencompeqother (ipp,0);
01459
01460
01461 for (j=0;j<nne;j++){
01462
01463 if (Tp->gradcomp==1){
01464 if (nodes[nod[j]].gradient==NULL){
01465 nodes[nod[j]].alloc_grad (ncomp);
01466 }
01467 }
01468
01469 if (Tp->fluxcomp==1){
01470 if (nodes[nod[j]].flux==NULL){
01471 nodes[nod[j]].alloc_flux (ncomp);
01472 }
01473 }
01474
01475 if (Tp->othercomp==1){
01476 if (nodes[nod[j]].other==NULL){
01477 nodes[nod[j]].alloc_other (ncompo);
01478 }
01479 }
01480
01481 if (Tp->eqothercomp==1){
01482 if (nodes[nod[j]].eqother==NULL){
01483 nodes[nod[j]].alloc_eqother (ncompeqo);
01484 }
01485 }
01486 }
01487 destrv (nod);
01488 }
01489
01490
01491 if (Tp->nvs==1){
01492 for (i=0;i<nn;i++){
01493 nodes[i].alloc_nodval ();
01494 }
01495 }
01496
01497 if (Tp->pnvs==1){
01498 for (i=0;i<nn;i++){
01499 nodes[i].alloc_nodvalp ();
01500 }
01501 }
01502
01503 if (Tp->invs==1){
01504 for (i=0;i<nn;i++){
01505 nodes[i].alloc_nodvali ();
01506 }
01507 }
01508
01509 if (Tp->tdnvs==1){
01510 for (i=0;i<nn;i++){
01511 nodes[i].alloc_nodvalt ();
01512 }
01513 }
01514 }
01515
01516
01517
01518
01519
01520
01521
01522 void transtop::alloc_enodes ()
01523 {
01524 if (Gtt->nen>0){
01525 enodes = new endnodet [Gtt->nen];
01526 }
01527 }
01528
01529
01530
01531
01532
01533
01534
01535 void transtop::alloc_edges ()
01536 {
01537 if (Gtt->nged>0){
01538 edges = new edget [Gtt->nged];
01539 }
01540 }
01541
01542
01543
01544
01545
01546
01547 void transtop::enodes_init ()
01548 {
01549 long i;
01550
01551 for (i=0;i<Gtt->nen;i++){
01552 enodes[i].init (i);
01553 }
01554 }
01555
01556
01557
01558
01559
01560
01561 void transtop::edge_init ()
01562 {
01563 long i;
01564
01565 for (i=0;i<Gtt->nged;i++){
01566 edges[i].init (i);
01567 }
01568 }
01569
01570
01571
01572
01573
01574
01575
01576 void transtop::edge_init_edval ()
01577 {
01578 long i;
01579
01580 for (i=0;i<Gtt->nged;i++){
01581 edges[i].init_edval ();
01582 }
01583 }
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594 void transtop::compute_jumps (double *rhs)
01595 {
01596 long i,j,nm,*ncn1,*ncn2,*mcn;
01597
01598
01599 for (i=0;i<Gtt->nen;i++){
01600 enodes[i].compute_jump (i);
01601
01602
01603 nm=Gtt->endnodes[i].ndofn;
01604 ncn1 = new long [nm];
01605 ncn2 = new long [nm];
01606 mcn = new long [nm];
01607
01608
01609 Gtt->give_endnode_code_numbers (i,ncn1,ncn2,mcn);
01610
01611 for (j=0;j<nm;j++){
01612 if (mcn[j]>0)
01613 rhs[mcn[j]-1]+=enodes[i].jump[j];
01614 }
01615 delete [] mcn;
01616 delete [] ncn2;
01617 delete [] ncn1;
01618 }
01619
01620
01621 for (i=0;i<Gtt->nged;i++){
01622 edges[i].compute_jump (i);
01623
01624
01625 nm=Gtt->gedges[i].ndofnf;
01626 ncn1 = new long [nm];
01627 ncn2 = new long [nm];
01628 mcn = new long [nm];
01629
01630
01631 Gtt->give_edge_code_numbers (i,1,ncn1,ncn2,mcn);
01632
01633 for (j=0;j<nm;j++){
01634 if (mcn[j]>0)
01635 rhs[mcn[j]-1]+=edges[i].jumpfn[j];
01636 }
01637 delete [] mcn;
01638 delete [] ncn2;
01639 delete [] ncn1;
01640
01641
01642
01643 nm=Gtt->gedges[i].ndofnl;
01644 ncn1 = new long [nm];
01645 ncn2 = new long [nm];
01646 mcn = new long [nm];
01647
01648
01649 Gtt->give_edge_code_numbers (i,2,ncn1,ncn2,mcn);
01650
01651 for (j=0;j<nm;j++){
01652 if (mcn[j]>0)
01653 rhs[mcn[j]-1]+=edges[i].jumpln[j];
01654 }
01655 delete [] mcn;
01656 delete [] ncn2;
01657 delete [] ncn1;
01658
01659 }
01660
01661 }
01662
01663
01664
01665
01666
01667
01668
01669 void transtop::compute_resistance_factor (double *rhs)
01670 {
01671 long i,j,k,edid,fnid,lnid,nm,lcid,*ncn1,*ncn2,*mcn;
01672 double ncf,*normal,*flux;
01673
01674
01675 lcid=0;
01676
01677
01678 normal = new double [2];
01679 flux = new double [2];
01680
01681
01682
01683 compute_nodefluxes ();
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720 for (i=0;i<Gtt->nser;i++){
01721
01722 for (j=0;j<Gtt->nedser[i];j++){
01723
01724 edid = Gtt->edgelist[i][j];
01725
01726
01727 fnid = Gtt->gedges[edid].fn;
01728
01729
01730 Gtt->gedges[edid].give_norvect (normal);
01731
01732
01733
01734 nodes[fnid].give_flux (lcid,flux);
01735
01736
01737
01738 ncf = ss (normal,flux,2);
01739
01740
01741
01742 edges[edid].compute_node_jump (edid,1,lcid,ncf);
01743
01744
01745 nm=Gtt->gedges[edid].ndofnf;
01746 ncn1 = new long [nm];
01747 ncn2 = new long [nm];
01748 mcn = new long [nm];
01749
01750
01751
01752 Gtt->give_edge_code_numbers (edid,1,ncn1,ncn2,mcn);
01753
01754 for (k=0;k<nm;k++){
01755 if (mcn[k]>0)
01756 rhs[mcn[k]-1]+=edges[edid].jumpfn[k];
01757 }
01758 delete [] mcn;
01759 delete [] ncn2;
01760 delete [] ncn1;
01761
01762
01763
01764 if (j==Gtt->nedser[i]){
01765
01766
01767 lnid = Gtt->gedges[edid].ln;
01768
01769
01770
01771 nodes[lnid].give_flux (lcid,flux);
01772
01773
01774
01775 ncf = ss (normal,flux,2);
01776
01777
01778
01779
01780 edges[edid].compute_node_jump (edid,2,lcid,ncf);
01781
01782
01783 nm=Gtt->gedges[edid].ndofnf;
01784 ncn1 = new long [nm];
01785 ncn2 = new long [nm];
01786 mcn = new long [nm];
01787
01788
01789 Gtt->give_edge_code_numbers (edid,2,ncn1,ncn2,mcn);
01790
01791 for (k=0;k<nm;k++){
01792 if (mcn[k]>0)
01793 rhs[mcn[k]-1]+=edges[edid].jumpln[k];
01794 }
01795 delete [] mcn;
01796 delete [] ncn2;
01797 delete [] ncn1;
01798
01799 }
01800 }
01801 }
01802
01803 delete [] flux;
01804 }
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814 void transtop::initial_nodval ()
01815 {
01816 long i,j,k,ndofe,lcid,*cn;
01817 double *r;
01818
01819
01820 lcid=0;
01821
01822 for (i=0;i<ne;i++){
01823
01824 j=Gtt->leso[i];
01825 k=Gtt->gelements[i].auxinf;
01826
01827 if (j==1 && k==0){
01828
01829
01830 ndofe = give_ndofe (i);
01831
01832 r = new double [ndofe];
01833 cn = new long [ndofe];
01834
01835
01836 give_code_numbers (i,cn);
01837
01838 elemvalues (lcid,i,r,cn,ndofe);
01839
01840 elements[i].initnodvalues (r,ndofe);
01841
01842 Gtt->gelements[i].auxinf=1;
01843
01844 delete [] cn;
01845 delete [] r;
01846 }
01847 }
01848
01849 }
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861 void transtop::lhs_save (double *lhs,double *lhsi,double *tdlhs)
01862 {
01863 long i,j,ndofn;
01864 long *cn;
01865
01866 for (i=0;i<nn;i++){
01867
01868 ndofn = give_ndofn (i);
01869 cn = new long [ndofn];
01870 give_node_code_numbers (i,cn);
01871
01872 for (j=0;j<ndofn;j++){
01873 nodes[i].nodval[j]=0.0;
01874 nodes[i].nodvali[j]=0.0;
01875 nodes[i].nodvalt[j]=0.0;
01876 if (cn[j]>0){
01877 nodes[i].nodval[j]=lhs[cn[j]-1];
01878 nodes[i].nodvali[j]=lhsi[cn[j]-1];
01879 nodes[i].nodvalt[j]=tdlhs[cn[j]-1];
01880 }
01881 }
01882 delete [] cn;
01883 }
01884
01885 }
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897 void transtop::lhs_restore (double *lhs,double *lhsi,double *tdlhs)
01898 {
01899 long i,j,ndofn;
01900 long *cn;
01901
01902 for (i=0;i<nn;i++){
01903 ndofn = give_ndofn (i);
01904 cn = new long [ndofn];
01905 give_node_code_numbers (i,cn);
01906
01907 for (j=0;j<ndofn;j++){
01908 if (cn[j]>0){
01909 lhs[cn[j]-1]=nodes[i].nodval[j];
01910 lhsi[cn[j]-1]=nodes[i].nodvali[j];
01911 tdlhs[cn[j]-1]=nodes[i].nodvalt[j];
01912 }
01913 }
01914 delete [] cn;
01915 }
01916
01917 }
01918
01919
01920 long transtop::mesh_check(void)
01921 {
01922 long i,j,err=0,kk,jj;
01923
01924
01925 for(i=0;i<ne;i++){
01926 for(j=0;j<Gtt->gelements[i].nne;j++){
01927
01928 jj = Gtt->gelements[i].nodes[j];
01929 for(kk = 0; kk < Tp->ntm;kk++){
01930
01931 if((Gtt->gelements[i].tgf) < (Gtt->gnodes[jj].tgf[kk])){
01932
01933 fprintf (stdout,"\n\n Different element node(%ld) time function number=%ld than element(%ld) time function number=%ld",j+1,Gtt->gnodes[jj].tgf[kk]+1,i+1,Gtt->gelements[i].tgf+1);
01934 fprintf (stdout,"\n Node number=%ld, element number=%ld\n\n",Gtt->gelements[i].nodes[j]+1,i+1);
01935 err = 1;
01936
01937 return err;
01938 }
01939
01940 }
01941 }
01942 }
01943
01944 return err;
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
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017 void transtop::view_factors (FILE *out)
02018 {
02019 long i, j, k;
02020 long emitter, emitter_fn, emitter_ln;
02021 long absorber, absorber_fn, absorber_ln;
02022 double diagonal1, diagonal2, side1, side2;
02023 double emitter_size,emitter_fn_Xcord, emitter_fn_Ycord, emitter_ln_Xcord, emitter_ln_Ycord;
02024 double absorber_fn_Xcord, absorber_fn_Ycord, absorber_ln_Xcord, absorber_ln_Ycord;
02025
02026
02027
02028 view_factor = new double** [Gtt->nser];
02029 for (i=0;i<Gtt->nser;i++){
02030 view_factor[i] = new double* [Gtt->nedser[i]];
02031 for (j=0;j<Gtt->nedser[i];j++){
02032 view_factor[i][j] = new double [Gtt->nedser[i]];
02033 }
02034 }
02035
02036
02037 for (i=0;i<Gtt->nser; i++){
02038
02039 for (j=0; j < Gtt->nedser[i] ; j++){
02040
02041 emitter = Gtt->edgelist[i][j];
02042
02043 emitter_fn = Gtt->gedges[emitter].fn;
02044 emitter_ln = Gtt->gedges[emitter].ln;
02045
02046
02047 emitter_fn_Xcord = Gtt->gnodes[emitter_fn].x;
02048 emitter_fn_Ycord = Gtt->gnodes[emitter_fn].y;
02049 emitter_ln_Xcord = Gtt->gnodes[emitter_ln].x;
02050 emitter_ln_Ycord = Gtt->gnodes[emitter_ln].y;
02051
02052
02053 emitter_size = sqrt(pow(emitter_fn_Xcord-emitter_ln_Xcord,2)+pow(emitter_fn_Ycord-emitter_ln_Ycord,2));
02054
02055 Gtt->gedges[emitter].l = emitter_size;
02056
02057
02058 for (k=0; k < Gtt->nedser[i]; k++) {
02059 absorber = Gtt->edgelist[i][k];
02060
02061 absorber_fn = Gtt->gedges[absorber].fn;
02062 absorber_ln = Gtt->gedges[absorber].ln;
02063
02064
02065 absorber_fn_Xcord = Gtt->gnodes[absorber_fn].x;
02066 absorber_fn_Ycord = Gtt->gnodes[absorber_fn].y;
02067 absorber_ln_Xcord = Gtt->gnodes[absorber_ln].x;
02068 absorber_ln_Ycord = Gtt->gnodes[absorber_ln].y;
02069
02070
02071 if (j == k) {
02072 view_factor[i][j][k] = 0;
02073 } else {
02074 diagonal1 = sqrt(pow(emitter_fn_Xcord-absorber_fn_Xcord,2)+pow(emitter_fn_Ycord-absorber_fn_Ycord,2));
02075 diagonal2 = sqrt(pow(emitter_ln_Xcord-absorber_ln_Xcord,2)+pow(emitter_ln_Ycord-absorber_ln_Ycord,2));
02076 side1 = sqrt(pow(emitter_fn_Xcord-absorber_ln_Xcord,2)+pow(emitter_fn_Ycord-absorber_ln_Ycord,2));
02077 side2 = sqrt(pow(emitter_ln_Xcord-absorber_fn_Xcord,2)+pow(emitter_ln_Ycord-absorber_fn_Ycord,2));
02078
02079 view_factor[i][j][k] = (diagonal1 + diagonal2 - side1 - side2)/(2*emitter_size);
02080 }
02081 }
02082 }
02083 }
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098 }
02099
02100
02101
02102
02103
02104
02105
02106
02107 void transtop::mod_view_factors (FILE *out)
02108 {
02109 long i,j,k,n,edid,e1id,e2id,matid;
02110 double eps,zero;
02111 double *x,*y;
02112 mattypet mat1,mat2;
02113 x=NULL;
02114 y=NULL;
02115
02116 zero=1.0e-15;
02117
02118 mvf = new densemat [Gtt->nser];
02119 mvfjk = new densemat [Gtt->nser];
02120
02121
02122 for (i=0;i<Gtt->nser; i++){
02123 n=Gtt->nedser[i];
02124 mvf[i].alloc (n);
02125 mvfjk[i].alloc (n);
02126
02127
02128 for (j=0;j<Gtt->nedser[i];j++){
02129
02130
02131 for (k=0; k < Gtt->nedser[i]; k++) {
02132
02133 edid = Gtt->edgelist[i][k];
02134
02135 e1id = Gtt->gedges[edid].adjel[0];
02136 e2id = Gtt->gedges[edid].adjel[1];
02137
02138 mat1 = elements[e1id].tm[0];
02139 mat2 = elements[e2id].tm[0];
02140
02141 if (mat1 == radiationmater){
02142
02143 matid = elements[e2id].idm[0];
02144
02145 eps = Tm->give_extinction_coeff (mat2,matid);
02146 }else{
02147
02148 matid = elements[e1id].idm[0];
02149
02150 eps = Tm->give_extinction_coeff (mat1,matid);
02151 }
02152
02153
02154 if (j==k){
02155 mvf[i].a[j*n+k]=1.0/eps;
02156 mvfjk[i].a[j*n+k]=1.0/eps;
02157 }else{
02158 mvf[i].a[j*n+k]=(1.0-1.0/eps)*view_factor[i][j][k];
02159 mvfjk[i].a[j*n+k]=(1.0-1.0/eps)*view_factor[i][j][k];
02160 }
02161 }
02162 }
02163
02164
02165 mvf[i].lu (x,y,zero,2);
02166 }
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181 }
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191 void transtop::edge_temperature ()
02192 {
02193 long i,j,edid,ndofn,fn,ln,lcid;
02194 double tfn,tln,taver;
02195 double *nval;
02196
02197
02198
02199 lcid=0;
02200
02201
02202 for (i=0;i<Gtt->nser;i++){
02203
02204 for (j=0;j<Gtt->nedser[i];j++){
02205
02206 edid = Gtt->edgelist[i][j];
02207
02208 fn = Gtt->gedges[edid].fn;
02209
02210 ln = Gtt->gedges[edid].ln;
02211
02212 ndofn = give_ndofn (fn);
02213
02214 nval = new double [ndofn];
02215
02216
02217 nodalval (lcid,nval,fn);
02218
02219 tfn = nval[0];
02220
02221 nodalval (lcid,nval,ln);
02222
02223 tln = nval[0];
02224
02225
02226 taver = (tfn+tln)/2.0;
02227
02228 edges[edid].store_edval (&taver);
02229
02230 delete [] nval;
02231 }
02232 }
02233 }
02234
02235
02236
02237
02238
02239
02240
02241
02242
02243 void transtop::heat_fluxes (double *rhs,FILE *out)
02244 {
02245 long i,j,edid,dof,fn,ln;
02246 double zero,l,h;
02247 double *x,*y,*z;
02248
02249 zero=1.0e-15;
02250
02251
02252 for (i=0;i<Gtt->nser;i++){
02253
02254 x = new double [Gtt->nedser[i]];
02255
02256 y = new double [Gtt->nedser[i]];
02257 z = new double [Gtt->nedser[i]];
02258
02259
02260 t4t4 (i,y);
02261
02262 mvf[i].lu (x,y,zero,3);
02263
02264 mvfjk[i].mxv_dm (x,z);
02265
02266 for (j=0;j<Gtt->nedser[i];j++){
02267
02268 edid = Gtt->edgelist[i][j];
02269
02270 l=Gtt->gedges[edid].l;
02271
02272 h=x[j]*l/2.0;
02273
02274 fn = Gtt->gedges[edid].fn;
02275
02276 ln = Gtt->gedges[edid].ln;
02277
02278 dof = give_dof (fn,0)-1;
02279 if (dof>-1){
02280 rhs[dof]-=h;
02281 }
02282 dof = give_dof (ln,0)-1;
02283 if (dof>-1){
02284 rhs[dof]-=h;
02285 }
02286
02287
02288 }
02289
02290 delete [] x;
02291 delete [] y;
02292 }
02293 }
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303 void transtop::t4t4 (long edserid,double *y)
02304 {
02305 long i,j,edid;
02306 double stefan_boltzmann,ti,tj;
02307
02308
02309 stefan_boltzmann=5.67e-8;
02310
02311
02312 for (i=0;i<Gtt->nedser[edserid];i++){
02313
02314 edid = Gtt->edgelist[edserid][i];
02315
02316 edges[edid].give_edval (&ti);
02317
02318 y[i]=0.0;
02319
02320
02321 for (j=0;j<Gtt->nedser[edserid];j++){
02322
02323 edid = Gtt->edgelist[edserid][j];
02324
02325 edges[edid].give_edval (&tj);
02326
02327 y[i]+=view_factor[edserid][i][j]*stefan_boltzmann*(ti*ti*ti*ti - tj*tj*tj*tj);
02328 }
02329 }
02330 }
02331