00001 #include <stdlib.h>
00002 #include "globmatt.h"
00003 #include "globalt.h"
00004 #include "lhsrhst.h"
00005 #include "constrel.h"
00006 #include "onemedium.h"
00007 #include "twomedia.h"
00008 #include "threemedia.h"
00009
00010
00011
00012
00013
00014 void conductivity_matrix (long lcid)
00015 {
00016 if (Kmat==NULL) Kmat = new gmatrix;
00017 Kmat->ts=(storagetype)(*(int*)&Tp->tstorkm);
00018 Kmat->initiate (Gtt,Ndoft,Mesprt);
00019
00020 Kmat->setval (Tp->ssle);
00021
00022 long i,ndofe,*cn;
00023 matrix lm;
00024
00025 for (i=0;i<Tt->ne;i++){
00026 ndofe = Gtt->give_ndofe (i);
00027
00028 allocm (ndofe,ndofe,lm);
00029 conductmat (i,lcid,lm);
00030
00031 cn = new long [ndofe];
00032 Gtt->give_code_numbers (i,cn);
00033
00034 Kmat->localize (lm,cn,i);
00035
00036 destrm (lm);
00037 delete [] cn;
00038 }
00039
00040 Kmat->prepmat (Gtt,0.0,Mesprt);
00041
00042 }
00043
00044
00045
00046
00047
00048
00049 void capacity_matrix (long lcid)
00050 {
00051 if (Cmat==NULL) Cmat = new gmatrix;
00052 Cmat->ts=(storagetype)(*(int*)&Tp->tstorcm);
00053 Cmat->initiate (Gtt,Ndoft,Mesprt);
00054
00055 Cmat->setval (Tp->ssle);
00056
00057 long i,ndofe,*cn;
00058 matrix lm;
00059
00060 for (i=0;i<Tt->ne;i++){
00061 ndofe = Gtt->give_ndofe (i);
00062
00063 allocm (ndofe,ndofe,lm);
00064
00065 capacmat (i,lcid,lm);
00066
00067 cn = new long [ndofe];
00068 Gtt->give_code_numbers (i,cn);
00069
00070 Cmat->localize (lm,cn,i);
00071
00072 destrm (lm);
00073 delete [] cn;
00074 }
00075
00076 Cmat->prepmat (Gtt,0.0,Mesprt);
00077
00078 }
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089 void residuum (double *r,double *p,double *v,double dt,long n,long lcid)
00090 {
00091 capacity_matrix (0);
00092 conductivity_matrix (0);
00093 trfel_right_hand_side (0,Lsrst->rhs,Ndoft);
00094
00095 Cmat->gmxv (p,v);
00096 cmulv (1.0/dt,v,n);
00097 Kmat->gmxv (Lsrst->lhs,r);
00098 addv (r,v,n);
00099 subv (r,Lsrst->rhs,n);
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 void nodalvalues (long lcid,double *r,long *cn,long ndofe)
00114 {
00115 long i,ii;
00116
00117 if ((Tp->tprob==nonstationary_problem) || (Tp->tprob==nonlinear_nonstationary_problem)){
00118 for (i=0;i<ndofe;i++){
00119 ii=cn[i];
00120 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00121 if (ii==0) r[i]=0.0;
00122 if (ii>0){
00123 r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 }
00145 }
00146 }
00147 else{
00148 for (i=0;i<ndofe;i++){
00149 ii=cn[i];
00150 if (ii<0) r[i]=Tb->lc[lcid].pv[0-ii-1].getval();
00151 if (ii==0) r[i]=0.0;
00152 if (ii>0) r[i]=Lsrst->lhs[ii-1];
00153 }
00154 }
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 }
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 void nodval (long lcid,double *r, long idn)
00190 {
00191 long i,ii;
00192
00193 if ((Tp->tprob==nonstationary_problem) || (Tp->tprob==nonlinear_nonstationary_problem)){
00194 for (i=0;i<Tt->give_ndofn (idn);i++){
00195 ii=Tt->give_dof(idn,i);
00196 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00197 if (ii==0) r[i]=0.0;
00198 if (ii>0) if (ii>0){
00199 r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218 }
00219 }
00220 }
00221 else{
00222 for (i=0;i<Tt->give_ndofn (idn);i++){
00223 ii=Tt->give_dof(idn,i);
00224 if (ii<0) r[i]=Tb->lc[lcid].pv[0-ii-1].getval();
00225 if (ii==0) r[i]=0.0;
00226 if (ii>0) r[i]=Lsrst->lhs[ii-1];
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 void nodalderivatives (double *r,long *cn,long ndofe)
00265 {
00266 long i,ii;
00267
00268 if ((Tp->tprob==nonstationary_problem) || (Tp->tprob==nonlinear_nonstationary_problem)){
00269 for (i=0;i<ndofe;i++){
00270 ii=cn[i];
00271 if (ii<0) r[i]=0.0;
00272 if (ii==0) r[i]=0.0;
00273 if (ii>0) r[i]=Lsrst->tdlhs[ii-1];
00274 }
00275 }
00276 }
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 void prescvalues (double *r,long *cn,long ndofe)
00288 {
00289 long i,ii;
00290
00291 for (i=0;i<ndofe;i++){
00292 ii=cn[i];
00293 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00294 if (ii==0) r[i]=0.0;
00295 if (ii>0) r[i]=Lsrst->lhsi[ii-1];
00296 }
00297 }
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308 void initialvalues (double *r,long *cn,long ndofe)
00309 {
00310 long i,ii;
00311
00312 for (i=0;i<ndofe;i++){
00313 ii=cn[i];
00314 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getipv();
00315 if (ii==0) r[i]=0.0;
00316 if (ii>0) r[i]=Lsrst->lhsi[ii-1];
00317 }
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 void conductmat (long eid,long lcid,matrix &km)
00329 {
00330 elemtypet te;
00331
00332 te=Tt->give_elem_type (eid);
00333
00334 switch (te){
00335 case barlint:{
00336 Lbt->res_conductivity_matrix (eid,lcid,km);
00337 break;
00338 }
00339 case barlintax:{
00340 Lbat->res_conductivity_matrix (eid,lcid,km);
00341 break;
00342 }
00343 case barquadt:{
00344 Qbt->res_conductivity_matrix (eid,lcid,km);
00345 break;
00346 }
00347 case barquadtax:{
00348 Qbat->res_conductivity_matrix (eid,lcid,km);
00349 break;
00350 }
00351 case trlint:{
00352 Ltt->res_conductivity_matrix (eid,lcid,km);
00353 break;
00354 }
00355 case trlaxisym:{
00356 Ltat->res_conductivity_matrix (eid,lcid,km);
00357 break;
00358 }
00359 case quadlint:{
00360 Lqt->res_conductivity_matrix (eid,lcid,km);
00361 break;
00362 }
00363 case quadquadt:{
00364 Qqt->res_conductivity_matrix (eid,lcid,km);
00365 break;
00366 }
00367 case quadquadtax:{
00368 Qqat->res_conductivity_matrix (eid,lcid,km);
00369 break;
00370 }
00371 case quadlaxisym:{
00372 Lqat->res_conductivity_matrix (eid,lcid,km);
00373 break;
00374 }
00375 case lineartett:{
00376 Ltett->res_conductivity_matrix (eid,lcid,km);
00377 break;
00378 }
00379 case linearhext:{
00380 Lht->res_conductivity_matrix (eid,lcid,km);
00381 break;
00382 }
00383 case quadratichext:{
00384 Qht->res_conductivity_matrix (eid,lcid,km);
00385 break;
00386 }
00387 default:{
00388 fprintf (stderr,"\n\n unknown element type is required in function");
00389 fprintf (stderr,"\n conductmat (file %s, line %d).\n",__FILE__,__LINE__);
00390 }
00391 }
00392
00393 }
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 void capacmat (long eid,long lcid,matrix &cm)
00404 {
00405 elemtypet te;
00406
00407 te=Tt->give_elem_type (eid);
00408
00409 switch (te){
00410 case barlint:{
00411 Lbt->res_capacity_matrix (eid,cm);
00412 break;
00413 }
00414 case barlintax:{
00415 Lbat->res_capacity_matrix (eid,cm);
00416 break;
00417 }
00418 case barquadt:{
00419 Qbt->res_capacity_matrix (eid,cm);
00420 break;
00421 }
00422 case barquadtax:{
00423 Qbat->res_capacity_matrix (eid,cm);
00424 break;
00425 }
00426 case trlint:{
00427 Ltt->res_capacity_matrix (eid,cm);
00428 break;
00429 }
00430 case trlaxisym:{
00431 Ltat->res_capacity_matrix (eid,cm);
00432 break;
00433 }
00434 case quadlint:{
00435 Lqt->res_capacity_matrix (eid,cm);
00436 break;
00437 }
00438 case quadquadt:{
00439 Qqt->res_capacity_matrix (eid,cm);
00440 break;
00441 }
00442 case quadquadtax:{
00443 Qqat->res_capacity_matrix (eid,cm);
00444 break;
00445 }
00446 case quadlaxisym:{
00447 Lqat->res_capacity_matrix (eid,cm);
00448 break;
00449 }
00450 case lineartett:{
00451 Ltett->res_capacity_matrix (eid,cm);
00452 break;
00453 }
00454 case linearhext:{
00455 Lht->res_capacity_matrix (eid,cm);
00456 break;
00457 }
00458 case quadratichext:{
00459 Qht->res_capacity_matrix (eid,cm);
00460 break;
00461 }
00462 default:{
00463 fprintf (stderr,"\n\n unknown element type is required in function");
00464 fprintf (stderr,"\n capacmat (file %s, line %d).\n",__FILE__,__LINE__);
00465 }
00466 }
00467 }
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 void internal_fluxes (double *intflux,long n)
00478 {
00479 long i,ne,ndofe,*cn;
00480 elemtypet te;
00481 vector lif;
00482
00483 ne = Tt->ne;
00484 nullv (intflux,n);
00485
00486 for (i=0;i<ne;i++){
00487 te = Tt->give_elem_type (i);
00488 ndofe = Gtt->give_ndofe (i);
00489 allocv (ndofe,lif);
00490 fillv (0.0,lif);
00491
00492 cn = new long [ndofe];
00493 Gtt->give_code_numbers (i,cn);
00494
00495 switch (te){
00496 case barlint:{
00497 Lbt->res_internal_fluxes (i,lif);
00498 break;
00499 }
00500 case barlintax:{
00501 Lbat->res_internal_fluxes (i,lif);
00502 break;
00503 }
00504 case barquadt:{
00505 Qbt->res_internal_fluxes (i,lif);
00506 break;
00507 }
00508 case barquadtax:{
00509 Qbat->res_internal_fluxes (i,lif);
00510 break;
00511 }
00512 case trlint:{
00513 Ltt->res_internal_fluxes (i,lif);
00514 break;
00515 }
00516 case trlaxisym:{
00517 Ltat->res_internal_fluxes (i,lif);
00518 break;
00519 }
00520 case quadlint:{
00521 Lqt->res_internal_fluxes (i,lif);
00522 break;
00523 }
00524 case quadquadt:{
00525 Qqt->res_internal_fluxes (i,lif);
00526 break;
00527 }
00528 case quadquadtax:{
00529 Qqat->res_internal_fluxes (i,lif);
00530 break;
00531 }
00532 case quadlaxisym:{
00533 Lqat->res_internal_fluxes (i,lif);
00534 break;
00535 }
00536 case lineartett:{
00537 Ltett->res_internal_fluxes (i,lif);
00538 break;
00539 }
00540 case linearhext:{
00541 Lht->res_internal_fluxes (i,lif);
00542 break;
00543 }
00544 case quadratichext:{
00545 Qht->res_internal_fluxes (i,lif);
00546 break;
00547 }
00548 default:{
00549 fprintf (stderr,"\n\n unknown element type is required in function internal_fluxes (file %s, line %d).\n",__FILE__,__LINE__);
00550 abort ();
00551 }
00552 }
00553
00554 locglob (intflux,lif.a,cn,ndofe);
00555 destrv (lif);
00556 delete [] cn;
00557
00558 }
00559
00560 trfel_bound_flux (0,intflux,n);
00561
00562 }
00563
00564
00565
00566
00567
00568
00569 void approximation ()
00570 {
00571 long i;
00572
00573 switch (Tp->tmatt){
00574 case onemedium:{
00575 for (i=0;i<Tt->ne;i++){
00576 intpointvalues (i);
00577 intpointgradients (i);
00578 }
00579 break;
00580 }
00581 case twomediacoup:{
00582 for (i=0;i<Tt->ne;i++){
00583 intpointvalues (i);
00584 intpointgradients (i);
00585 }
00586 break;
00587 }
00588 case threemediacoup:{
00589 for (i=0;i<Tt->ne;i++){
00590 intpointvalues (i);
00591 intpointgradients (i);
00592 }
00593 break;
00594 }
00595 default:{
00596 fprintf (stderr,"\n\n unknown type of element in function approximation (file %s, line %d).\n",__FILE__,__LINE__);
00597 }
00598 }
00599
00600 }
00601
00602 void intpointvalues (long eid)
00603 {
00604 elemtypet te;
00605
00606 te=Tt->give_elem_type (eid);
00607
00608 switch (te){
00609 case barlint:{
00610 Lbt->intpointval (eid);
00611 break;
00612 }
00613 case barlintax:{
00614 Lbat->intpointval (eid);
00615 break;
00616 }
00617 case barquadt:{
00618 Qbt->intpointval (eid);
00619 break;
00620 }
00621 case barquadtax:{
00622 Qbat->intpointval (eid);
00623 break;
00624 }
00625 case trlint:{
00626 Ltt->intpointval (eid);
00627 break;
00628 }
00629 case trlaxisym:{
00630 Ltat->intpointval (eid);
00631 break;
00632 }
00633 case quadlint:{
00634 Lqt->intpointval (eid);
00635 break;
00636 }
00637 case quadquadt:{
00638 Qqt->intpointval (eid);
00639 break;
00640 }
00641 case quadquadtax:{
00642 Qqat->intpointval (eid);
00643 break;
00644 }
00645 case quadlaxisym:{
00646 Lqat->intpointval (eid);
00647 break;
00648 }
00649 case lineartett:{
00650 Ltett->intpointval (eid);
00651 break;
00652 }
00653 case linearhext:{
00654 Lht->intpointval (eid);
00655 break;
00656 }
00657 case quadratichext:{
00658 Qht->intpointval (eid);
00659 break;
00660 }
00661 default:{
00662 fprintf (stderr,"\n\n unknown element type is required in function intpointvalues (file %s, line %d).\n",__FILE__,__LINE__);
00663 }
00664 }
00665
00666 }
00667
00668 void intpointgradients (long eid)
00669 {
00670 elemtypet te;
00671
00672 te=Tt->give_elem_type (eid);
00673
00674 switch (te){
00675 case barlint:{
00676 Lbt->intpointgrad (eid);
00677 break;
00678 }
00679 case barlintax:{
00680 Lbat->intpointgrad (eid);
00681 break;
00682 }
00683 case barquadt:{
00684 Qbt->intpointgrad (eid);
00685 break;
00686 }
00687 case barquadtax:{
00688 Qbat->intpointgrad (eid);
00689 break;
00690 }
00691 case trlint:{
00692 Ltt->intpointgrad (eid);
00693 break;
00694 }
00695 case trlaxisym:{
00696 Ltat->intpointgrad (eid);
00697 break;
00698 }
00699 case quadlint:{
00700 Lqt->intpointgrad (eid);
00701 break;
00702 }
00703 case quadquadt:{
00704 Qqt->intpointgrad (eid);
00705 break;
00706 }
00707 case quadquadtax:{
00708 Qqat->intpointgrad (eid);
00709 break;
00710 }
00711 case quadlaxisym:{
00712 Lqat->intpointgrad (eid);
00713 break;
00714 }
00715 case lineartett:{
00716 Ltett->intpointgrad (eid);
00717 break;
00718 }
00719 case linearhext:{
00720 Lht->intpointgrad (eid);
00721 break;
00722 }
00723 case quadratichext:{
00724 Qht->intpointgrad (eid);
00725 break;
00726 }
00727 default:{
00728 fprintf (stderr,"\n\n unknown element type is required in function intpointvalues (file %s, line %d).\n",__FILE__,__LINE__);
00729 }
00730 }
00731
00732 }
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742 void assemble_init (double *rhs)
00743 {
00744 long i,ndofe,*cn,lcid=0,ii,kk;
00745 vector r,f;
00746 matrix km;
00747 long ne = Tt->ne;
00748
00749 for (i=0;i<ne;i++){
00750 ndofe=Tt->give_ndofe(i);
00751 allocm (ndofe,ndofe,km);
00752
00753 conductmat (i,lcid,km);
00754
00755 allocv (ndofe,r);
00756 allocv (ndofe,f);
00757 cn = new long [ndofe];
00758
00759 Tt->give_code_numbers (i,cn);
00760
00761 prescvalues (r.a,cn,ndofe);
00762
00763
00764 mxv (km,r,f);
00765
00766 cmulv (-1.0,f);
00767
00768 locglob (rhs,f.a,cn,ndofe);
00769
00770 destrm (km); destrv (f); destrv (r);
00771 delete [] cn;
00772 }
00773 }
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783 void trfel_right_hand_side (long lcid,double *rhs,long n)
00784 {
00785 double *av;
00786 av = new double [n];
00787 nullv (av,n);
00788
00789 nullv (rhs,n);
00790
00791 assemble_init (rhs);
00792
00793 switch (Tp->tmatt){
00794 case onemedium:{
00795 Tb->lc[lcid].assemble (lcid,av+lcid*n);
00796 addv (rhs,av,n);
00797 break;
00798 }
00799 case twomediacoup:{
00800 Tb->lc[0].assemble (0,av);
00801 addv (rhs,av,n);
00802 nullv (av,n);
00803 Tb->lc[1].assemble (1,av);
00804 addv (rhs,av,n);
00805 break;
00806 }
00807 case threemediacoup:{
00808 Tb->lc[0].assemble (0,av);
00809 addv (rhs,av,n);
00810 nullv (av,n);
00811 Tb->lc[1].assemble (1,av);
00812 addv (rhs,av,n);
00813 nullv (av,n);
00814 Tb->lc[2].assemble (2,av);
00815 addv (rhs,av,n);
00816 break;
00817 }
00818 default:{
00819 fprintf (stderr,"\n\n unknown transported matter is required in function");
00820 fprintf (stderr,"\n right_hand_side (file %s, line %d).\n",__FILE__,__LINE__);
00821 }
00822 }
00823
00824 delete [] av;
00825 }
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836 void trfel_bound_flux (long lcid,double *iflux,long n)
00837 {
00838 double *av;
00839 av = new double [n];
00840
00841 switch (Tp->tmatt){
00842 case onemedium:{
00843 Tb->lc[lcid].assemble_flux (lcid,av+lcid*n,n);
00844 addv (iflux,av,n);
00845 break;
00846 }
00847 case twomediacoup:{
00848 Tb->lc[0].assemble_flux (0,av,n);
00849 addv (iflux,av,n);
00850 nullv (av,n);
00851 Tb->lc[1].assemble_flux (1,av,n);
00852 addv (iflux,av,n);
00853 break;
00854 }
00855 case threemediacoup:{
00856 Tb->lc[0].assemble_flux (0,av,n);
00857 addv (iflux,av,n);
00858 nullv (av,n);
00859 Tb->lc[1].assemble_flux (1,av,n);
00860 addv (iflux,av,n);
00861 nullv (av,n);
00862 Tb->lc[2].assemble_flux (2,av,n);
00863 addv (iflux,av,n);
00864 break;
00865 }
00866 default:{
00867 fprintf (stderr,"\n\n unknown transported matter is required in function");
00868 fprintf (stderr,"\n right_hand_side (file %s, line %d).\n",__FILE__,__LINE__);
00869 }
00870 }
00871
00872 delete [] av;
00873 }
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886 void compute_req_valt (long lcid)
00887 {
00888
00889
00890
00891
00892
00893
00894
00895 if (Tp->othercomp==1)
00896 Tm->compute_nodeothers (lcid);
00897 }
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909 void source_vector (long lcid,long eid,vector &nodval,vector &lv)
00910 {
00911 elemtypet te;
00912
00913
00914 te=Tt->elements[eid].te;
00915
00916 switch (te){
00917 case barlint:{
00918 Lbt->res_quantity_source_vector (lv,nodval,lcid,eid);
00919 break;
00920 }
00921 case barlintax:{
00922 Lbat->res_quantity_source_vector (lv,nodval,lcid,eid);
00923 break;
00924 }
00925 case trlint:{
00926 Ltt->res_quantity_source_vector (lv,nodval,lcid,eid);
00927 break;
00928 }
00929 case trlaxisym:{
00930 Ltat->res_quantity_source_vector (lv,nodval,lcid,eid);
00931 break;
00932 }
00933 case quadlint:{
00934 Lqt->res_quantity_source_vector (lv,nodval,lcid,eid);
00935 break;
00936 }
00937 case quadlaxisym:{
00938 Lqat->res_quantity_source_vector (lv,nodval,lcid,eid);
00939 break;
00940 }
00941 case quadquadt:{
00942 Qqt->res_quantity_source_vector (lv,nodval,lcid,eid);
00943 break;
00944 }
00945 case quadquadtax:{
00946 Qqat->res_quantity_source_vector (lv,nodval,lcid,eid);
00947 break;
00948 }
00949 case lineartett:{
00950 Ltett->res_quantity_source_vector (lv,nodval,lcid,eid);
00951 break;
00952 }
00953 case linearhext:{
00954 Lht->res_quantity_source_vector (lv,nodval,lcid,eid);
00955 break;
00956 }
00957 case quadratichext:{
00958 Qht->res_quantity_source_vector (lv,nodval,lcid,eid);
00959 break;
00960 }
00961 default:{
00962 fprintf (stderr,"\n\n unknown element type is required in function source_vector (file %s, line %d",__FILE__,__LINE__);
00963 }
00964 }
00965
00966 }
00967
00968
00969 void give_nodal_humid (double *gv,long mnt)
00970 {
00971 long i,cn,compother;
00972 double *r;
00973
00974 r = new double [Tp->ntm];
00975
00976 switch (Tp->tmatt){
00977 case onemedium:{
00978 med1 m1;
00979
00980 for (i=0;i<mnt;i++){
00981 cn=Tt->give_dof (i,0);
00982 if (cn>0) gv[i] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00983 if (cn<0) gv[i] = Tb->lc[0].pv[0-cn-1].getval();
00984 }
00985 break;
00986
00987 gv[i] = m1.compute_othervalues(1,0,r);
00988 if(gv[i] >= 1.0)
00989 gv[i] = 1.0;
00990 }
00991 case twomediacoup:{
00992 med2 m2;
00993
00994 for (i=0;i<mnt;i++){
00995 cn=Tt->give_dof (i,0);
00996 if (cn>0) r[0] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00997 if (cn==0) r[0] = 0.0;
00998 if (cn<0) r[0] = Tb->lc[0].pv[0-cn-1].getval();
00999
01000 cn=Tt->give_dof (i,1);
01001 if (cn>0) r[1] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01002 if (cn==0) r[1] = 0.0;
01003 if (cn<0) r[1] = Tb->lc[0].pv[0-cn-1].getval();
01004
01005 compother = 1;
01006 if ((Tm->ip[0].tm) == kunzel)
01007 compother = 0;
01008
01009 gv[i] = m2.compute_othervalues(compother,0,r);
01010 if(gv[i] >= 1.0)
01011 gv[i] = 1.0;
01012 }
01013 break;
01014 }
01015 case threemediacoup:{
01016 med3 m3;
01017
01018 for (i=0;i<mnt;i++){
01019 cn=Tt->give_dof (i,0);
01020 if (cn>0) r[0] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01021 if (cn==0) r[0] = 0.0;
01022 if (cn<0) r[0] = Tb->lc[0].pv[0-cn-1].getval();
01023
01024 cn=Tt->give_dof (i,1);
01025 if (cn>0) r[1] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01026 if (cn==0) r[1] = 0.0;
01027 if (cn<0) r[1] = Tb->lc[0].pv[0-cn-1].getval();
01028
01029 cn=Tt->give_dof (i,2);
01030 if (cn>0) r[2] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
01031 if (cn==0) r[2] = 0.0;
01032 if (cn<0) r[2] = Tb->lc[0].pv[0-cn-1].getval();
01033
01034 gv[i] = m3.compute_othervalues(3,0,r);
01035 if(gv[i] >= 1.0)
01036 gv[i] = 1.0;
01037 }
01038 break;
01039 }
01040 default:{
01041 fprintf (stderr,"\n\n unknown type of transported matter in function (file %s, line %d).\n",__FILE__,__LINE__);
01042 }
01043 }
01044
01045 delete [] r;
01046 }
01047
01048
01049
01050
01051
01052
01053
01054
01055 void solution_correction ()
01056 {
01057 long i,j,k,cn,ne,nne,ndofn,lcid,ipp;
01058 ivector nodes;
01059 vector nv;
01060
01061
01062 ne = Tt->ne;
01063
01064
01065 lcid=0;
01066
01067
01068 for (i=0;i<ne;i++){
01069
01070 nne = Tt->give_nne (i);
01071
01072
01073 allocv (nne,nodes);
01074
01075
01076 Tt->give_elemnodes (i,nodes);
01077
01078
01079 ipp=Tt->elements[i].ipp[0][0];
01080
01081
01082 for (j=0;j<nne;j++){
01083
01084 ndofn = Tt->give_ndofn (nodes[j]);
01085
01086 allocv (ndofn,nv);
01087
01088 nodval (lcid,nv.a,nodes[j]);
01089
01090 Tm->values_correction (nv,ipp);
01091
01092
01093 for (k=0;k<ndofn;k++){
01094 cn = Tt->give_dof (nodes[j],k);
01095 if (cn>0){
01096 Lsrst->lhs[cn-1]=nv[k]-Lsrst->lhsi[cn-1];
01097 }
01098 }
01099
01100
01101 destrv (nv);
01102 }
01103
01104 destrv (nodes);
01105 }
01106
01107 }