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) r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00123 }
00124 }
00125 else{
00126 for (i=0;i<ndofe;i++){
00127 ii=cn[i];
00128 if (ii<0) r[i]=Tb->lc[lcid].pv[0-ii-1].getval();
00129 if (ii==0) r[i]=0.0;
00130 if (ii>0) r[i]=Lsrst->lhs[ii-1];
00131 }
00132 }
00133
00134 }
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 void nodval (long lcid,double *r, long idn)
00147 {
00148 long i,ii;
00149
00150 if ((Tp->tprob==nonstationary_problem) || (Tp->tprob==nonlinear_nonstationary_problem)){
00151 for (i=0;i<Tt->give_ndofn (idn);i++){
00152 ii=Tt->give_dof(idn,i);
00153 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00154 if (ii==0) r[i]=0.0;
00155 if (ii>0) r[i]=Lsrst->lhsi[ii-1]+Lsrst->lhs[ii-1];
00156 }
00157 }
00158 else{
00159 for (i=0;i<Tt->give_ndofn (idn);i++){
00160 ii=Tt->give_dof(idn,i);
00161 if (ii<0) r[i]=Tb->lc[lcid].pv[0-ii-1].getval();
00162 if (ii==0) r[i]=0.0;
00163 if (ii>0) r[i]=Lsrst->lhs[ii-1];
00164 }
00165 }
00166
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 void nodalderivatives (double *r,long *cn,long ndofe)
00180 {
00181 long i,ii;
00182
00183 if ((Tp->tprob==nonstationary_problem) || (Tp->tprob==nonlinear_nonstationary_problem)){
00184 for (i=0;i<ndofe;i++){
00185 ii=cn[i];
00186 if (ii<0) r[i]=0.0;
00187 if (ii==0) r[i]=0.0;
00188 if (ii>0) r[i]=Lsrst->tdlhs[ii-1];
00189 }
00190 }
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 void prescvalues (double *r,long *cn,long ndofe)
00203 {
00204 long i,ii;
00205
00206 for (i=0;i<ndofe;i++){
00207 ii=cn[i];
00208 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getval();
00209 if (ii==0) r[i]=0.0;
00210 if (ii>0) r[i]=Lsrst->lhsi[ii-1];
00211 }
00212 }
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223 void initialvalues (double *r,long *cn,long ndofe)
00224 {
00225 long i,ii;
00226
00227 for (i=0;i<ndofe;i++){
00228 ii=cn[i];
00229 if (ii<0) r[i]=Tb->lc[0].pv[0-ii-1].getipv();
00230 if (ii==0) r[i]=0.0;
00231 if (ii>0) r[i]=Lsrst->lhsi[ii-1];
00232 }
00233 }
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 void conductmat (long eid,long lcid,matrix &km)
00244 {
00245 elemtypet te;
00246
00247 te=Tt->give_elem_type (eid);
00248
00249 switch (te){
00250 case barlint:{
00251 Lbt->res_conductivity_matrix (eid,lcid,km);
00252 break;
00253 }
00254 case barlintax:{
00255 Lbat->res_conductivity_matrix (eid,lcid,km);
00256 break;
00257 }
00258 case barquadt:{
00259 Qbt->res_conductivity_matrix (eid,lcid,km);
00260 break;
00261 }
00262 case barquadtax:{
00263 Qbat->res_conductivity_matrix (eid,lcid,km);
00264 break;
00265 }
00266 case trlint:{
00267 Ltt->res_conductivity_matrix (eid,lcid,km);
00268 break;
00269 }
00270 case trlaxisym:{
00271 Ltat->res_conductivity_matrix (eid,lcid,km);
00272 break;
00273 }
00274 case quadlint:{
00275 Lqt->res_conductivity_matrix (eid,lcid,km);
00276 break;
00277 }
00278 case quadquadt:{
00279 Qqt->res_conductivity_matrix (eid,lcid,km);
00280 break;
00281 }
00282 case quadquadtax:{
00283 Qqat->res_conductivity_matrix (eid,lcid,km);
00284 break;
00285 }
00286 case quadlaxisym:{
00287 Lqat->res_conductivity_matrix (eid,lcid,km);
00288 break;
00289 }
00290 case lineartett:{
00291 Ltett->res_conductivity_matrix (eid,lcid,km);
00292 break;
00293 }
00294 case linearhext:{
00295 Lht->res_conductivity_matrix (eid,lcid,km);
00296 break;
00297 }
00298 case quadratichext:{
00299 Qht->res_conductivity_matrix (eid,lcid,km);
00300 break;
00301 }
00302 default:{
00303 fprintf (stderr,"\n\n unknown element type is required in function");
00304 fprintf (stderr,"\n conductmat (file %s, line %d).\n",__FILE__,__LINE__);
00305 }
00306 }
00307
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318 void capacmat (long eid,long lcid,matrix &cm)
00319 {
00320 elemtypet te;
00321
00322 te=Tt->give_elem_type (eid);
00323
00324 switch (te){
00325 case barlint:{
00326 Lbt->res_capacity_matrix (eid,cm);
00327 break;
00328 }
00329 case barlintax:{
00330 Lbat->res_capacity_matrix (eid,cm);
00331 break;
00332 }
00333 case barquadt:{
00334 Qbt->res_capacity_matrix (eid,cm);
00335 break;
00336 }
00337 case barquadtax:{
00338 Qbat->res_capacity_matrix (eid,cm);
00339 break;
00340 }
00341 case trlint:{
00342 Ltt->res_capacity_matrix (eid,cm);
00343 break;
00344 }
00345 case trlaxisym:{
00346 Ltat->res_capacity_matrix (eid,cm);
00347 break;
00348 }
00349 case quadlint:{
00350 Lqt->res_capacity_matrix (eid,cm);
00351 break;
00352 }
00353 case quadquadt:{
00354 Qqt->res_capacity_matrix (eid,cm);
00355 break;
00356 }
00357 case quadquadtax:{
00358 Qqat->res_capacity_matrix (eid,cm);
00359 break;
00360 }
00361 case quadlaxisym:{
00362 Lqat->res_capacity_matrix (eid,cm);
00363 break;
00364 }
00365 case lineartett:{
00366 Ltett->res_capacity_matrix (eid,cm);
00367 break;
00368 }
00369 case linearhext:{
00370 Lht->res_capacity_matrix (eid,cm);
00371 break;
00372 }
00373 case quadratichext:{
00374 Qht->res_capacity_matrix (eid,cm);
00375 break;
00376 }
00377 default:{
00378 fprintf (stderr,"\n\n unknown element type is required in function");
00379 fprintf (stderr,"\n capacmat (file %s, line %d).\n",__FILE__,__LINE__);
00380 }
00381 }
00382 }
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392 void internal_fluxes (double *intflux,long n)
00393 {
00394 long i,ne,ndofe,*cn;
00395 elemtypet te;
00396 vector lif;
00397
00398 ne = Tt->ne;
00399 nullv (intflux,n);
00400
00401 for (i=0;i<ne;i++){
00402 te = Tt->give_elem_type (i);
00403 ndofe = Gtt->give_ndofe (i);
00404 allocv (ndofe,lif);
00405 fillv (0.0,lif);
00406
00407 cn = new long [ndofe];
00408 Gtt->give_code_numbers (i,cn);
00409
00410 switch (te){
00411 case barlint:{
00412 Lbt->res_internal_fluxes (i,lif);
00413 break;
00414 }
00415 case barlintax:{
00416 Lbat->res_internal_fluxes (i,lif);
00417 break;
00418 }
00419 case barquadt:{
00420 Qbt->res_internal_fluxes (i,lif);
00421 break;
00422 }
00423 case barquadtax:{
00424 Qbat->res_internal_fluxes (i,lif);
00425 break;
00426 }
00427 case trlint:{
00428 Ltt->res_internal_fluxes (i,lif);
00429 break;
00430 }
00431 case trlaxisym:{
00432 Ltat->res_internal_fluxes (i,lif);
00433 break;
00434 }
00435 case quadlint:{
00436 Lqt->res_internal_fluxes (i,lif);
00437 break;
00438 }
00439 case quadquadt:{
00440 Qqt->res_internal_fluxes (i,lif);
00441 break;
00442 }
00443 case quadquadtax:{
00444 Qqat->res_internal_fluxes (i,lif);
00445 break;
00446 }
00447 case quadlaxisym:{
00448 Lqat->res_internal_fluxes (i,lif);
00449 break;
00450 }
00451 case lineartett:{
00452 Ltett->res_internal_fluxes (i,lif);
00453 break;
00454 }
00455 case linearhext:{
00456 Lht->res_internal_fluxes (i,lif);
00457 break;
00458 }
00459 case quadratichext:{
00460 Qht->res_internal_fluxes (i,lif);
00461 break;
00462 }
00463 default:{
00464 fprintf (stderr,"\n\n unknown element type is required in function internal_fluxes (file %s, line %d).\n",__FILE__,__LINE__);
00465 abort ();
00466 }
00467 }
00468
00469 locglob (intflux,lif.a,cn,ndofe);
00470 destrv (lif);
00471 delete [] cn;
00472
00473 }
00474
00475 trfel_bound_flux (0,intflux,n);
00476
00477 }
00478
00479
00480
00481
00482
00483
00484 void approximation ()
00485 {
00486 long i;
00487
00488 switch (Tp->tmatt){
00489 case onemedium:{
00490 for (i=0;i<Tt->ne;i++){
00491 intpointvalues (i);
00492 intpointgradients (i);
00493 }
00494 break;
00495 }
00496 case twomediacoup:{
00497 for (i=0;i<Tt->ne;i++){
00498 intpointvalues (i);
00499 intpointgradients (i);
00500 }
00501 break;
00502 }
00503 case threemediacoup:{
00504 for (i=0;i<Tt->ne;i++){
00505 intpointvalues (i);
00506 intpointgradients (i);
00507 }
00508 break;
00509 }
00510 default:{
00511 fprintf (stderr,"\n\n unknown type of element in function approximation (file %s, line %d).\n",__FILE__,__LINE__);
00512 }
00513 }
00514
00515 }
00516
00517 void intpointvalues (long eid)
00518 {
00519 elemtypet te;
00520
00521 te=Tt->give_elem_type (eid);
00522
00523 switch (te){
00524 case barlint:{
00525 Lbt->intpointval (eid);
00526 break;
00527 }
00528 case barlintax:{
00529 Lbat->intpointval (eid);
00530 break;
00531 }
00532 case barquadt:{
00533 Qbt->intpointval (eid);
00534 break;
00535 }
00536 case barquadtax:{
00537 Qbat->intpointval (eid);
00538 break;
00539 }
00540 case trlint:{
00541 Ltt->intpointval (eid);
00542 break;
00543 }
00544 case trlaxisym:{
00545 Ltat->intpointval (eid);
00546 break;
00547 }
00548 case quadlint:{
00549 Lqt->intpointval (eid);
00550 break;
00551 }
00552 case quadquadt:{
00553 Qqt->intpointval (eid);
00554 break;
00555 }
00556 case quadquadtax:{
00557 Qqat->intpointval (eid);
00558 break;
00559 }
00560 case quadlaxisym:{
00561 Lqat->intpointval (eid);
00562 break;
00563 }
00564 case lineartett:{
00565 Ltett->intpointval (eid);
00566 break;
00567 }
00568 case linearhext:{
00569 Lht->intpointval (eid);
00570 break;
00571 }
00572 case quadratichext:{
00573 Qht->intpointval (eid);
00574 break;
00575 }
00576 default:{
00577 fprintf (stderr,"\n\n unknown element type is required in function intpointvalues (file %s, line %d).\n",__FILE__,__LINE__);
00578 }
00579 }
00580
00581 }
00582
00583 void intpointgradients (long eid)
00584 {
00585 elemtypet te;
00586
00587 te=Tt->give_elem_type (eid);
00588
00589 switch (te){
00590 case barlint:{
00591 Lbt->intpointgrad (eid);
00592 break;
00593 }
00594 case barlintax:{
00595 Lbat->intpointgrad (eid);
00596 break;
00597 }
00598 case barquadt:{
00599 Qbt->intpointgrad (eid);
00600 break;
00601 }
00602 case barquadtax:{
00603 Qbat->intpointgrad (eid);
00604 break;
00605 }
00606 case trlint:{
00607 Ltt->intpointgrad (eid);
00608 break;
00609 }
00610 case trlaxisym:{
00611 Ltat->intpointgrad (eid);
00612 break;
00613 }
00614 case quadlint:{
00615 Lqt->intpointgrad (eid);
00616 break;
00617 }
00618 case quadquadt:{
00619 Qqt->intpointgrad (eid);
00620 break;
00621 }
00622 case quadquadtax:{
00623 Qqat->intpointgrad (eid);
00624 break;
00625 }
00626 case quadlaxisym:{
00627 Lqat->intpointgrad (eid);
00628 break;
00629 }
00630 case lineartett:{
00631 Ltett->intpointgrad (eid);
00632 break;
00633 }
00634 case linearhext:{
00635 Lht->intpointgrad (eid);
00636 break;
00637 }
00638 case quadratichext:{
00639 Qht->intpointgrad (eid);
00640 break;
00641 }
00642 default:{
00643 fprintf (stderr,"\n\n unknown element type is required in function intpointvalues (file %s, line %d).\n",__FILE__,__LINE__);
00644 }
00645 }
00646
00647 }
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657 void assemble_init (double *rhs)
00658 {
00659 long i,ndofe,*cn,lcid=0,ii,kk;
00660 vector r,f;
00661 matrix km;
00662 long ne = Tt->ne;
00663
00664 for (i=0;i<ne;i++){
00665 ndofe=Tt->give_ndofe(i);
00666 allocm (ndofe,ndofe,km);
00667
00668 conductmat (i,lcid,km);
00669
00670 allocv (ndofe,r);
00671 allocv (ndofe,f);
00672 cn = new long [ndofe];
00673
00674 Tt->give_code_numbers (i,cn);
00675
00676 prescvalues (r.a,cn,ndofe);
00677
00678
00679 mxv (km,r,f);
00680
00681 cmulv (-1.0,f);
00682
00683 locglob (rhs,f.a,cn,ndofe);
00684
00685 destrm (km); destrv (f); destrv (r);
00686 delete [] cn;
00687 }
00688 }
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698 void trfel_right_hand_side (long lcid,double *rhs,long n)
00699 {
00700 double *av;
00701 av = new double [n];
00702 nullv (av,n);
00703
00704 nullv (rhs,n);
00705
00706 assemble_init (rhs);
00707
00708 switch (Tp->tmatt){
00709 case onemedium:{
00710 Tb->lc[lcid].assemble (lcid,av+lcid*n);
00711 addv (rhs,av,n);
00712 break;
00713 }
00714 case twomediacoup:{
00715 Tb->lc[0].assemble (0,av);
00716 addv (rhs,av,n);
00717 nullv (av,n);
00718 Tb->lc[1].assemble (1,av);
00719 addv (rhs,av,n);
00720 break;
00721 }
00722 case threemediacoup:{
00723 Tb->lc[0].assemble (0,av);
00724 addv (rhs,av,n);
00725 nullv (av,n);
00726 Tb->lc[1].assemble (1,av);
00727 addv (rhs,av,n);
00728 nullv (av,n);
00729 Tb->lc[2].assemble (2,av);
00730 addv (rhs,av,n);
00731 break;
00732 }
00733 default:{
00734 fprintf (stderr,"\n\n unknown transported matter is required in function");
00735 fprintf (stderr,"\n right_hand_side (file %s, line %d).\n",__FILE__,__LINE__);
00736 }
00737 }
00738
00739 delete [] av;
00740 }
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751 void trfel_bound_flux (long lcid,double *iflux,long n)
00752 {
00753 double *av;
00754 av = new double [n];
00755
00756 switch (Tp->tmatt){
00757 case onemedium:{
00758 Tb->lc[lcid].assemble_flux (lcid,av+lcid*n,n);
00759 addv (iflux,av,n);
00760 break;
00761 }
00762 case twomediacoup:{
00763 Tb->lc[0].assemble_flux (0,av,n);
00764 addv (iflux,av,n);
00765 nullv (av,n);
00766 Tb->lc[1].assemble_flux (1,av,n);
00767 addv (iflux,av,n);
00768 break;
00769 }
00770 case threemediacoup:{
00771 Tb->lc[0].assemble_flux (0,av,n);
00772 addv (iflux,av,n);
00773 nullv (av,n);
00774 Tb->lc[1].assemble_flux (1,av,n);
00775 addv (iflux,av,n);
00776 nullv (av,n);
00777 Tb->lc[2].assemble_flux (2,av,n);
00778 addv (iflux,av,n);
00779 break;
00780 }
00781 default:{
00782 fprintf (stderr,"\n\n unknown transported matter is required in function");
00783 fprintf (stderr,"\n right_hand_side (file %s, line %d).\n",__FILE__,__LINE__);
00784 }
00785 }
00786
00787 delete [] av;
00788 }
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801 void compute_req_valt (long lcid)
00802 {
00803
00804
00805
00806
00807
00808
00809
00810 if (Tp->othercomp==1)
00811 Tm->compute_nodeothers (lcid);
00812 }
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824 void source_vector (long lcid,long eid,vector &nodval,vector &lv)
00825 {
00826 elemtypet te;
00827
00828
00829 te=Tt->elements[eid].te;
00830
00831 switch (te){
00832 case barlint:{
00833 Lbt->res_quantity_source_vector (lv,nodval,lcid,eid);
00834 break;
00835 }
00836 case barlintax:{
00837 Lbat->res_quantity_source_vector (lv,nodval,lcid,eid);
00838 break;
00839 }
00840 case trlint:{
00841 Ltt->res_quantity_source_vector (lv,nodval,lcid,eid);
00842 break;
00843 }
00844 case trlaxisym:{
00845 Ltat->res_quantity_source_vector (lv,nodval,lcid,eid);
00846 break;
00847 }
00848 case quadlint:{
00849 Lqt->res_quantity_source_vector (lv,nodval,lcid,eid);
00850 break;
00851 }
00852 case quadlaxisym:{
00853 Lqat->res_quantity_source_vector (lv,nodval,lcid,eid);
00854 break;
00855 }
00856 case quadquadt:{
00857 Qqt->res_quantity_source_vector (lv,nodval,lcid,eid);
00858 break;
00859 }
00860 case quadquadtax:{
00861 Qqat->res_quantity_source_vector (lv,nodval,lcid,eid);
00862 break;
00863 }
00864 case lineartett:{
00865 Ltett->res_quantity_source_vector (lv,nodval,lcid,eid);
00866 break;
00867 }
00868 case linearhext:{
00869 Lht->res_quantity_source_vector (lv,nodval,lcid,eid);
00870 break;
00871 }
00872 case quadratichext:{
00873 Qht->res_quantity_source_vector (lv,nodval,lcid,eid);
00874 break;
00875 }
00876 default:{
00877 fprintf (stderr,"\n\n unknown element type is required in function source_vector (file %s, line %d",__FILE__,__LINE__);
00878 }
00879 }
00880
00881 }
00882
00883
00884 void give_nodal_humid (double *gv,long mnt)
00885 {
00886 long i,cn,compother;
00887 double *r;
00888
00889 r = new double [Tp->ntm];
00890
00891 switch (Tp->tmatt){
00892 case onemedium:{
00893 med1 m1;
00894
00895 for (i=0;i<mnt;i++){
00896 cn=Tt->give_dof (i,0);
00897 if (cn>0) gv[i] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00898 if (cn<0) gv[i] = Tb->lc[0].pv[0-cn-1].getval();
00899 }
00900 break;
00901
00902 gv[i] = m1.compute_othervalues(1,0,r);
00903 if(gv[i] >= 1.0)
00904 gv[i] = 1.0;
00905 }
00906 case twomediacoup:{
00907 med2 m2;
00908
00909 for (i=0;i<mnt;i++){
00910 cn=Tt->give_dof (i,0);
00911 if (cn>0) r[0] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00912 if (cn==0) r[0] = 0.0;
00913 if (cn<0) r[0] = Tb->lc[0].pv[0-cn-1].getval();
00914
00915 cn=Tt->give_dof (i,1);
00916 if (cn>0) r[1] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00917 if (cn==0) r[1] = 0.0;
00918 if (cn<0) r[1] = Tb->lc[0].pv[0-cn-1].getval();
00919
00920 compother = 1;
00921 if ((Tm->ip[0].tm) == kunzel)
00922 compother = 0;
00923
00924 gv[i] = m2.compute_othervalues(compother,0,r);
00925 if(gv[i] >= 1.0)
00926 gv[i] = 1.0;
00927 }
00928 break;
00929 }
00930 case threemediacoup:{
00931 med3 m3;
00932
00933 for (i=0;i<mnt;i++){
00934 cn=Tt->give_dof (i,0);
00935 if (cn>0) r[0] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00936 if (cn==0) r[0] = 0.0;
00937 if (cn<0) r[0] = Tb->lc[0].pv[0-cn-1].getval();
00938
00939 cn=Tt->give_dof (i,1);
00940 if (cn>0) r[1] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00941 if (cn==0) r[1] = 0.0;
00942 if (cn<0) r[1] = Tb->lc[0].pv[0-cn-1].getval();
00943
00944 cn=Tt->give_dof (i,2);
00945 if (cn>0) r[2] = Lsrst->lhsi[cn-1]+Lsrst->lhs[cn-1];
00946 if (cn==0) r[2] = 0.0;
00947 if (cn<0) r[2] = Tb->lc[0].pv[0-cn-1].getval();
00948
00949 gv[i] = m3.compute_othervalues(3,0,r);
00950 if(gv[i] >= 1.0)
00951 gv[i] = 1.0;
00952 }
00953 break;
00954 }
00955 default:{
00956 fprintf (stderr,"\n\n unknown type of transported matter in function (file %s, line %d).\n",__FILE__,__LINE__);
00957 }
00958 }
00959
00960 delete [] r;
00961 }
00962
00963
00964
00965
00966
00967
00968
00969 void solution_correction ()
00970 {
00971 long i,j,k,cn,ne,nne,ndofn,lcid,ipp;
00972 ivector nodes;
00973 vector nv;
00974
00975
00976 ne = Tt->ne;
00977
00978
00979 lcid=0;
00980
00981
00982 for (i=0;i<ne;i++){
00983
00984 nne = Tt->give_nne (i);
00985
00986
00987 allocv (nne,nodes);
00988
00989
00990 Tt->give_elemnodes (i,nodes);
00991
00992
00993 ipp=Tt->elements[i].ipp[0][0];
00994
00995
00996 for (j=0;j<nne;j++){
00997
00998 ndofn = Tt->give_ndofn (nodes[j]);
00999
01000 allocv (ndofn,nv);
01001
01002 nodval (lcid,nv.a,nodes[j]);
01003
01004 Tm->values_correction (nv,ipp);
01005
01006
01007 for (k=0;k<ndofn;k++){
01008 cn = Tt->give_dof (nodes[j],k);
01009 if (cn>0){
01010 Lsrst->lhs[cn-1]=nv[k]-Lsrst->lhsi[cn-1];
01011 }
01012 }
01013
01014
01015 destrv (nv);
01016 }
01017
01018 destrv (nodes);
01019 }
01020
01021 }