00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "kunzel2.h"
00029 #include "globalt.h"
00030
00031
00032 kunmat2::kunmat2()
00033 {
00034 rho_m = 0.0;
00035 rho_w = 1000.0;
00036 }
00037
00038 kunmat2::~kunmat2()
00039 {
00040
00041
00042 }
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 void kunmat2::matcond (matrix &d,long ri,long ci,long ipp)
00057 {
00058 long n;
00059 n = d.n;
00060
00061 switch (n){
00062 case 1:{
00063 matcond1d (d,ri,ci,ipp);
00064 break;
00065 }
00066 case 2:{
00067 matcond2d (d,ri,ci,ipp);
00068 break;
00069 }
00070 case 3:{
00071 matcond3d (d,ri,ci,ipp);
00072 break;
00073 }
00074 default:{
00075 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00076 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00077 }
00078 }
00079 }
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090 void kunmat2::matcond1d (matrix &d,long ri,long ci,long ipp)
00091 {
00092 double k;
00093 double x1,x2;
00094 k = 0.0;
00095
00096 x1 = Tm->ip[ipp].av[0];
00097 x2 = Tm->ip[ipp].av[1];
00098
00099 if((ri == 0) && (ci == 0))
00100 k = tokJ1(x1,x2,ipp);
00101 if((ri == 0) && (ci == 1))
00102 k = tokJ2(x1,x2,ipp);
00103 if((ri == 1) && (ci == 0))
00104 k = tokJ3(x1,x2,ipp);
00105 if((ri == 1) && (ci == 1))
00106 k = tokJ4(x1,x2,ipp);
00107
00108 d[0][0] = k;
00109 }
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 void kunmat2::matcond2d (matrix &d,long ri,long ci,long ipp)
00120 {
00121 double k;
00122 double x1,x2;
00123 k = 0.0;
00124
00125 x1 = Tm->ip[ipp].av[0];
00126 x2 = Tm->ip[ipp].av[1];
00127
00128 if((ri == 0) && (ci == 0))
00129 k = tokJ1(x1,x2,ipp);
00130 if((ri == 0) && (ci == 1))
00131 k = tokJ2(x1,x2,ipp);
00132 if((ri == 1) && (ci == 0))
00133 k = tokJ3(x1,x2,ipp);
00134 if((ri == 1) && (ci == 1))
00135 k = tokJ4(x1,x2,ipp);
00136
00137 fillm(0.0,d);
00138
00139 d[0][0] = k; d[0][1] = 0.0;
00140 d[1][0] = 0.0; d[1][1] = k;
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 void kunmat2::matcond3d (matrix &d,long ri,long ci,long ipp)
00153 {
00154 double k;
00155 double x1,x2;
00156 k = 0.0;
00157
00158 x1 = Tm->ip[ipp].av[0];
00159 x2 = Tm->ip[ipp].av[1];
00160
00161 if((ri == 0) && (ci == 0))
00162 k = tokJ1(x1,x2,ipp);
00163 if((ri == 0) && (ci == 1))
00164 k = tokJ2(x1,x2,ipp);
00165 if((ri == 1) && (ci == 0))
00166 k = tokJ3(x1,x2,ipp);
00167 if((ri == 1) && (ci == 1))
00168 k = tokJ4(x1,x2,ipp);
00169
00170 fillm(0.0,d);
00171
00172 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00173 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00174 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00175 }
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 void kunmat2::matcap (double &c,long ri,long ci,long ipp)
00187 {
00188 double x1,x2;
00189 c = 0.0;
00190
00191 x1 = Tm->ip[ipp].av[0];
00192 x2 = Tm->ip[ipp].av[1];
00193
00194 if((ri == 0) && (ci == 0))
00195
00196 c = Tm->ip[ipp].eqother[1];
00197 if((ri == 0) && (ci == 1))
00198 c = 0.0;
00199 if((ri == 1) && (ci == 0))
00200 c = 0.0;
00201 if((ri == 1) && (ci == 1))
00202 c = DerivaceHustotyEntalpiePodleTeploty(x1,x2, ipp);
00203 }
00204
00205
00206 void kunmat2::values_correction (vector &nv, long ipp)
00207 {
00208
00209 relhum_check(nv[0],nv[1],ipp);
00210 }
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222 void kunmat2::relhum_check(double &x1,double x2,long ipp)
00223 {
00224 if (x1 >= 1.0)
00225 x1 = 1.0;
00226 if (x1 <= 0.0)
00227 x1 = 0.0;
00228
00229
00230 Tm->ip[ipp].av[0] = x1;
00231 }
00232
00233
00234
00235
00236
00237
00238
00239
00240 double kunmat2::DerivaceTlakuNasycenychParNaTeplote(double x1, double x2)
00241 {
00242
00243 return ((4042.9 * exp(23.5771-4042.9/(x2 - 37.58)))/((x2-37.58)*(x2-37.58)));
00244 }
00245
00246
00247
00248
00249
00250
00251
00252
00253 double kunmat2::PermeabilitaVodniPary(double x1, double x2, long ipp)
00254 {
00255 double P, Rv, Pa;
00256
00257 double da, dp, mi, smc, moist1;
00258
00259 if (x1 > 1) x1 = 1;
00260 if (x2 < 0) x1 = 0;
00261
00262
00263
00264 moist1 = Tm->ip[ipp].eqother[0];
00265 smc = Tm->ip[ipp].eqother[2];
00266 CorD(4,kd,smc,moist1,mi,a2,a3);
00267
00268 Rv = 461.5;
00269 Pa = 101325;
00270 P = 101325;
00271
00272 da = (2.306e-5 * Pa)/(Rv * x2 * P)*pow((x2/273.15),1.81);
00273 dp = da/mi;
00274 return (dp);
00275 }
00276
00277
00278
00279
00280
00281
00282 double kunmat2::tokJ2(double x1, double x2, long ipp)
00283 {
00284
00285 double dp, dpsdt;
00286
00287
00288 if (x1 > 1) x1 = 1;
00289 if (x1 < 0) x1 = 0;
00290
00291
00292 dp = DerivaceTlakuNasycenychParNaTeplote(x1,x2);
00293 dpsdt = PermeabilitaVodniPary(x1,x2, ipp);
00294
00295 J2 = dp * x1 * dpsdt;
00296
00297
00298
00299 return (J2);
00300
00301 }
00302
00303
00304
00305
00306
00307
00308
00309
00310 double kunmat2::tokJ4(double x1, double x2, long ipp)
00311 {
00312 double lv1, lambdaakt, dpsdt, dp, moist1,smc;
00313
00314
00315 if (x1 > 1) x1 = 1;
00316 if (x1 < 0) x1 = 0;
00317
00318
00319 moist1 = Tm->ip[ipp].eqother[0];
00320 smc = Tm->ip[ipp].eqother[2];
00321
00322 CorD(10,kd,smc,moist1,lambdaakt,a2,a3);
00323
00324 dpsdt = DerivaceTlakuNasycenychParNaTeplote(x1, x2);
00325 dp = PermeabilitaVodniPary(x1, x2, ipp);
00326
00327 lv1 = LatentHeatofEvaporationOfWater (x1, x2);
00328
00329 J4 = lambdaakt + lv1 * dp * x1 * dpsdt;
00330
00331 return (J4);
00332
00333 }
00334
00335
00336
00337
00338
00339
00340
00341
00342 double kunmat2::tokJ3(double x1, double x2, long ipp)
00343 {
00344 double lv1;
00345 double dp, ps;
00346
00347 if (x1 > 1) x1 = 1;
00348 if (x1 < 0) x1 = 0;
00349
00350 dp = PermeabilitaVodniPary(x1, x2, ipp);
00351 ps = TlakNasycenychVodnichParNaTeplote (x1, x2);
00352 lv1 = LatentHeatofEvaporationOfWater (x1, x2);
00353
00354 J3 = lv1 * dp * ps;
00355
00356 return (J3);
00357 }
00358
00359
00360 double kunmat2::LatentHeatofEvaporationOfWater(double x1, double x2)
00361 {
00362 return (2.5008e6)*pow((273.15/x2),(0.167+x2*3.67e-4));
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377 double kunmat2::TlakNasycenychVodnichParNaTeplote(double x1, double x2)
00378 {
00379 return (exp(23.5771 - 4042.9/(x2 - 37.58)));
00380 }
00381
00382
00383
00384
00385
00386
00387
00388
00389 double kunmat2::DerivaceHustotyEntalpiePodleTeploty(double x1, double x2, long ipp)
00390 {
00391 double rom, c, moist1, smc;
00392 if (x1 > 1) x1 = 1;
00393 if (x1 < 0) x1 = 0;
00394
00395 moist1 = Tm->ip[ipp].eqother[0];
00396
00397
00398 smc = Tm->ip[ipp].eqother[2];
00399 CorD(2,kd,0,0,rom,a2,a3);
00400
00401 rho_m = rom;
00402
00403 CorD(9,kd,smc,moist1,c,a2,a3);
00404
00405 return (rom * c);
00406 }
00407
00408
00409
00410
00411
00412
00413
00414 double kunmat2::HygroscopicMoisture(double x1, double x2)
00415 {
00416
00417 double mhmc, mhrh;
00418 CorD(6,kd,0,x1,mhmc,mhrh,a3);
00419
00420 return (mhmc/(1-sqrt(1-mhrh)));
00421
00422 }
00423
00424
00425
00426
00427
00428
00429
00430
00431 double kunmat2::DerivativeOfTheSorptionIsotherm(double x1, double x2)
00432 {
00433
00434 double mhv;
00435 if (x1 > 1) x1 = 1;
00436 if (x1 < 0) x1 = 0;
00437
00438 mhv = HygroscopicMoisture (x1,x2);
00439
00440 return (1000 * mhv/(2*sqrt(1-x1)));
00441
00442 }
00443
00444
00445
00446
00447
00448
00449
00450 double kunmat2::DerivativeOfTheRetentionCurve(double x1, double x2)
00451 {
00452
00453 double smc, mhmc, mhrh;
00454 CorD(7,kd,0,x1,smc,a2,a3);
00455 CorD(6,kd,0,x2,mhmc,mhrh,a3);
00456
00457
00458 return (1000 * (smc - mhmc)/(1 - mhrh));
00459 }
00460
00461
00462
00463
00464
00465
00466
00467 double kunmat2::DerivativeOfTheMoistureRetentionCharacteristik(double x1,double x2, long ipp)
00468 {
00469
00470 double drfdrh;
00471 drfdrh = 1e20;
00472
00473
00474 return (drfdrh);
00475
00476 }
00477
00478
00479
00480
00481
00482
00483
00484 double kunmat2::tokJ1(double x1, double x2, long ipp)
00485 {
00486
00487 double kapa1, mi, dp, dmretc, ps, drh, moist1, smc;
00488
00489
00490
00491 if (x1 < 0) x1 = 0;
00492 if (x1 > 1) x1 = 1;
00493
00494
00495 moist1 = Tm->ip[ipp].eqother[0];
00496 smc = Tm->ip[ipp].eqother[2];
00497
00498 CorD(4,kd,smc,moist1,mi,a2,a3);
00499
00500
00501
00502 kapa1 = Tm->ip[ipp].eqother[3];
00503
00504
00505 dp = PermeabilitaVodniPary(x1, x2, ipp);
00506 ps = TlakNasycenychVodnichParNaTeplote (x1, x2);
00507
00508 dmretc = Tm->ip[ipp].eqother[1];
00509 drh = kapa1 * dmretc;
00510 J1 = drh + dp * ps;
00511
00512 return (J1);
00513 }
00514
00515 void kunmat2::read(XFILE *in)
00516 {
00517 int i, j;
00518
00519 for (i = 2; i<19; i++)
00520 {
00521 xfscanf(in, "%d",&MatChar[i]);
00522 }
00523 xfscanf(in, "%d",&MatChar[19]);
00524 for (i = 2;i <= 19; i++)
00525 {
00526 switch (MatChar[i]){
00527 case 0:{
00528 break;
00529 }
00530 case 1:{
00531 xfscanf(in, "%lf",&MatConst[i]);
00532 break;
00533 }
00534 case 2:{
00535 if( MatChar[i-1] ==2)
00536 {
00537 xfscanf(in, "%lf %lf",&MatData[i][0][0],&MatData[i][1][0]);
00538 }
00539 else
00540 {
00541 xfscanf(in, "%lf %lf",&MatData[i][0][0],&MatData[i][1][0]);
00542 }
00543 switch (int(MatData[i][1][0])){
00544 case 2:{
00545 for (j=1; j<=MatData[i][0][0];j++)
00546 {
00547 xfscanf(in, "%lf %lf",&MatData[i][0][j],&MatData[i][1][j] );
00548 }
00549 break;
00550 }
00551 case 3:{
00552 for (j=1; j<=MatData[i][0][0];j++)
00553 {
00554 xfscanf(in, "%lf %lf %lf",&MatData[i][0][j],&MatData[i][1][j],&MatData[i][2][j] );
00555 }
00556 break;
00557 }
00558 default:{
00559 fprintf (stderr,"\n\n unknown definition of Material KUNZEL-DATA is required");
00560 fprintf (stderr,"\n in function read (file %s, line %d)\n",__FILE__,__LINE__);
00561 }
00562 }
00563 xfscanf(in, "%d",&MatData[i][0][0]);
00564 break;
00565 }
00566 case 3:{
00567 break;
00568 }
00569 case 30:{
00570 if( MatChar[i-1] ==1)
00571 {
00572 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00573 }
00574 else
00575 {
00576 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00577 }
00578 break;
00579 }
00580 case 31:{
00581 if( MatChar[i-1] ==1)
00582 {
00583 if (i == 6)
00584 {
00585 xfscanf(in, "%lf %lf %lf",&MatFunce[i][0],&MatFunce[i][1], &MatFunce[i][2]);
00586 }
00587 else
00588 {
00589 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00590 }
00591 }
00592 else
00593 {
00594 if (i == 6)
00595 {
00596 xfscanf(in, "%lf %lf %lf",&MatFunce[i][0],&MatFunce[i][1], &MatFunce[i][2]);
00597 }
00598 else
00599 {
00600 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00601 }
00602 }
00603 break;
00604 }
00605 case 32:{
00606 if( MatChar[i-1] ==1)
00607 {
00608 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00609 }
00610 else
00611 {
00612 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00613 }
00614 break;
00615 }
00616 case 40:{
00617
00618 xfscanf(in, "%lf %lf",&MatData[i][0][0],&Init[i]);
00619 xfscanf(in, "%lf %lf %lf %lf %lf %lf",&k1[i], &k2[i], &k3[i],&k4[i], &k5[i], &k6[i]);
00620 for (j=1; j<=MatData[i][0][0];j++)
00621 {
00622 xfscanf(in, "%lf %lf",&MatData[i][0][j],&MatData[i][1][j]);
00623 }
00624 xfscanf(in, "%lf",&MatData[i][0][0]);
00625 xfscanf(in, "%lf",&MatData[i][2][0]);
00626 for (j=1; j<=MatData[i][2][0];j++)
00627 {
00628 xfscanf(in, "%lf %lf",&MatData[i][2][j] ,&MatData[i][3][j] );
00629 }
00630 xfscanf(in, "%lf",&MatData[i][2][0]);
00631 break;
00632 }
00633 default:{
00634 fprintf (stderr,"\n\n unknown definition of Material KUNZEL is required");
00635 fprintf (stderr,"\n in function read (file %s, line %d)\n",__FILE__,__LINE__);
00636 }
00637 }
00638 }
00639 }
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767 double kunmat2::soptionizothermDerivation(double x1, double x2)
00768 {
00769 double a, b, n;
00770
00771 CorD(6,kd,0,x1,a,b,n);
00772
00773
00774 if (x1 > 1) x1 = 1;
00775 if (x1 < 0) x1 = 0;
00776
00777
00778 return (1000 * (a/(b*x1*n))*pow((1-(log(x1))/b),(-1-(1/n))));
00779
00780 }
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794 double kunmat2::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp)
00795 {
00796 long k;
00797 double new_trc,h,t;
00798 new_trc = 0.0;
00799
00800 k=Gtt->give_dof(nn,0);
00801 if (k>0) {h = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00802 if (k==0) {h = 0.0;}
00803 if (k<0) {h = Tb->lc[0].pv[0-k-1].getval();}
00804 k=Gtt->give_dof(nn,1);
00805 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00806 if (k==0) {t = 0.0;}
00807 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00808
00809 if((ri == 0) && (ci == 0))
00810 new_trc = get_transmission_transcoeff_hh(h,t,bc,ipp);
00811 if((ri == 0) && (ci == 1))
00812 new_trc = 0.0;
00813
00814 if((ri == 1) && (ci == 0))
00815 new_trc = 0.0;
00816 if((ri == 1) && (ci == 1))
00817 new_trc = get_transmission_transcoeff_tt(h,t,bc,ipp);
00818
00819 new_trc = new_trc*trc;
00820
00821 return (new_trc);
00822 }
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837 double kunmat2::transmission_nodval(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
00838 {
00839 long k;
00840 double new_nodval,h,t;
00841 new_nodval = 0.0;
00842
00843 k=Gtt->give_dof(nn,0);
00844 if (k>0) {h = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00845 if (k==0) {h = 0.0;}
00846 if (k<0) {h = Tb->lc[0].pv[0-k-1].getval();}
00847 k=Gtt->give_dof(nn,1);
00848 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00849 if (k==0) {t = 0.0;}
00850 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00851
00852 if((ri == 0) && (ci == 0))
00853 new_nodval = get_transmission_nodval_hh(nodval,h,t,bc,ipp);
00854 if((ri == 0) && (ci == 1))
00855 new_nodval = 0.0;
00856
00857 if((ri == 1) && (ci == 0))
00858 new_nodval = 0.0;
00859 if((ri == 1) && (ci == 1))
00860 new_nodval = get_transmission_nodval_tt(nodval,h,t,bc,ipp);
00861
00862 return (new_nodval);
00863 }
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879 double kunmat2::transmission_flux(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
00880 {
00881 long k;
00882 double flux,h,t;
00883 flux = 0.0;
00884
00885 k=Gtt->give_dof(nn,0);
00886 if (k>0) {h = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00887 if (k==0) {h = 0.0;}
00888 if (k<0) {h = Tb->lc[0].pv[0-k-1].getval();}
00889 k=Gtt->give_dof(nn,1);
00890 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00891 if (k==0) {t = 0.0;}
00892 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00893
00894 if((ri == 0) && (ci == 0))
00895 flux = get_transmission_flux_hh(nodval,h,t,bc,ipp);
00896 if((ri == 0) && (ci == 1))
00897 flux = 0.0;
00898
00899 if((ri == 1) && (ci == 0))
00900 flux = 0.0;
00901 if((ri == 1) && (ci == 1))
00902 flux = get_transmission_flux_tt(nodval,h,t,bc,ipp);
00903
00904 return (flux);
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915 double kunmat2::get_transmission_transcoeff_hh(double x1,double x2,long bc,long ipp)
00916 {
00917 double trc,pgws;
00918
00919 switch (bc){
00920 case 30:{
00921 pgws = TlakNasycenychVodnichParNaTeplote(x1,x2);
00922 trc = pgws;
00923 break;
00924 }
00925 case 31:{
00926 pgws = TlakNasycenychVodnichParNaTeplote(x1,x2);
00927 trc = pgws;
00928 break;
00929 }
00930 case 32:{
00931 trc = 0.0;
00932 break;
00933 }
00934 default:{
00935 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
00936 exit(0);
00937 }
00938 }
00939 return(trc);
00940 }
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950 double kunmat2::get_transmission_nodval_hh(double bv,double x1,double x2,long bc,long ipp)
00951 {
00952 double nodval,pgws, pgwspred;
00953
00954 switch (bc){
00955 case 30:{
00956
00957 pgws = TlakNasycenychVodnichParNaTeplote(x1,x2);
00958 pgwspred = TlakNasycenychVodnichParNaTeplote(x1,x2);
00959
00960 nodval = bv*pgwspred;
00961
00962 break;
00963 }
00964 case 31:{
00965 nodval = bv;
00966 break;
00967 }
00968 case 32:{
00969 pgws = TlakNasycenychVodnichParNaTeplote(x1,x2);
00970
00971 nodval = pgws*x1;
00972 bv = pgws*bv;
00973 nodval = bv - nodval;
00974
00975 break;
00976 }
00977 default:{
00978 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
00979 exit(0);
00980 }
00981 }
00982 return(nodval);
00983 }
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994 double kunmat2::get_transmission_flux_hh(double bv,double x1,double x2,long bc,long ipp)
00995 {
00996 double flux,pgws;
00997
00998 switch (bc){
00999 case 30:{
01000 pgws = TlakNasycenychVodnichParNaTeplote(x1, x2);
01001
01002 flux = pgws*x1;
01003 bv = pgws*bv;
01004 flux = bv - flux;
01005
01006 break;
01007 }
01008 case 31:{
01009 pgws = TlakNasycenychVodnichParNaTeplote(x1, x2);
01010 flux = pgws*x1;
01011
01012 flux = bv - flux;
01013 break;
01014 }
01015 case 32:{
01016 pgws = TlakNasycenychVodnichParNaTeplote(x1, x2);
01017
01018 flux = pgws*x1;
01019 bv = pgws*bv;
01020 flux = bv - flux;
01021
01022 break;
01023 }
01024 default:{
01025 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01026 exit(0);
01027 }
01028 }
01029 return(flux);
01030 }
01031
01032
01033
01034
01035
01036
01037
01038
01039 double kunmat2::get_transmission_transcoeff_tt(double x1,double x2,long bc,long ipp)
01040 {
01041 double trc;
01042
01043 switch (bc){
01044 case 30:{
01045 trc = 1.0;
01046 break;
01047 }
01048 default:{
01049 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01050 exit(0);
01051 }
01052 }
01053 return(trc);
01054 }
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064 double kunmat2::get_transmission_nodval_tt(double bv,double x1,double x2,long bc,long ipp)
01065 {
01066 double nodval;
01067
01068 switch (bc){
01069 case 30:{
01070 nodval = bv;
01071 break;
01072 }
01073 default:{
01074 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01075 exit(0);
01076 }
01077 }
01078 return(nodval);
01079 }
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089 double kunmat2::get_transmission_flux_tt(double bv,double x1,double x2,long bc,long ipp)
01090 {
01091 double flux;
01092
01093 switch (bc){
01094 case 30:{
01095 flux = (bv - x2);
01096 break;
01097 }
01098 default:{
01099 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01100 exit(0);
01101 }
01102 }
01103 return(flux);
01104 }
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115 double kunmat2::get_othervalue(long compother, double x1,double x2, long ipp)
01116 {
01117 double other;
01118 double t;
01119 t = Tp->time;
01120
01121
01122 switch (compother){
01123 case 0:{
01124 other = x1;
01125 break;
01126 }
01127 case 1:{
01128 other = x2;
01129 break;
01130 }
01131 case 2:{
01132
01133 other = Tm->ip[ipp].eqother[0];
01134
01135
01136
01137
01138
01139
01140
01141
01142 break;
01143 }
01144 case 3:{
01145 other = x1 * TlakNasycenychVodnichParNaTeplote (x1, x2);
01146 break;
01147 }
01148 case 4:{
01149 double w;
01150 w = Tm->ip[ipp].eqother[0];
01151 other = w*rho_w/rho_m;
01152 break;
01153 }
01154 case 5:other = Tm->ip[ipp].eqother[2];
01155 case 6:
01156 case 7:
01157 break;
01158
01159 default:{
01160 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
01161 }
01162 }
01163 return (other);
01164
01165 }
01166
01167
01168
01169
01170
01171
01172 void kunmat2::print_othervalue_name(FILE *out,long compother)
01173 {
01174 switch (compother){
01175 case 0:{
01176 fprintf (out,"Relative humidity () ");
01177 break;
01178 }
01179 case 1:{
01180 fprintf (out,"Temperature (K) ");
01181 break;
01182 }
01183 case 2:{
01184 fprintf (out,"Moisture content w (m3/m3) ");
01185 break;
01186 }
01187 case 3:{
01188 fprintf (out,"Water vapour pressure (Pa) ");
01189 break;
01190 }
01191 case 4:{
01192 fprintf (out,"Moisture content u (kg/kg) ");
01193 break;
01194 }
01195 default:{
01196 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
01197 }
01198 }
01199 }
01200
01201
01202 void kunmat2::give_data(double x1,double Mhmc, double Smc, double Mhrh, double & moistakt)
01203 {
01204 double hvezdicka;
01205 hvezdicka = Mhmc/(1-sqrt(1-Mhrh));
01206
01207 if (x1 < Mhrh) {
01208 moistakt = (1-sqrt(1-x1)) * hvezdicka ;
01209 }
01210 else {
01211 if (x1 > 1)
01212 {
01213 moistakt = Smc;
01214 }
01215 else
01216 {
01217 moistakt = Mhmc + (x1 - Mhrh )/(1-Mhrh )* (Smc - Mhmc );
01218 }
01219 }
01220
01221 }
01222
01223
01224 double kunmat2::kapa_exp(double a, double b, double x1w, double x2, long ipp)
01225 {
01226 if (x1w <= 0) x1w = 0;
01227 return (a * exp(b * x1w));
01228 }
01229
01230
01231 void kunmat2::CorD(int cislochar, int &kvyhl,double in, double x, double & y, double & z, double &z2)
01232 {
01233 int i;
01234 switch (MatChar[cislochar]){
01235 case 0:{
01236 kvyhl = 0;
01237 break;
01238 }
01239 case 1:{
01240 y = MatConst[cislochar];
01241 z = -10000;
01242 z2 = -10000;
01243 kvyhl = 1;
01244 break;
01245 }
01246 case 2:{
01247 int i = 1;
01248
01249 int psl, prad;
01250 prad = int( MatData[cislochar][0][0]);
01251 psl = int (MatData[cislochar][1][0]);
01252 if (x <=0)
01253 {
01254 x = 0;
01255 y = MatData[cislochar][1][1];
01256 if (psl == 3)
01257 {
01258 z = MatData[cislochar][2][1];
01259 }
01260
01261 }
01262 else {
01263 if (x > in)
01264 {
01265 x = in;
01266 y = MatData[cislochar][1][prad];
01267 if (psl == 3)
01268 {
01269 z = MatData[cislochar][2][prad];
01270 }
01271 }
01272 else
01273 {
01274 while (x > MatData[cislochar][0][i]){
01275 i++;
01276 }
01277
01278 y = MatData[cislochar][1][i-1] + ((x - MatData[cislochar][0][i-1])*(MatData[cislochar][1][i] - MatData[cislochar][1][i-1]))/(MatData[cislochar][0][i] - MatData[cislochar][0][i-1]);
01279 if (psl == 3)
01280 {
01281 z = MatData[cislochar][2][i-1] + ((x - MatData[cislochar][0][i-1])*(MatData[cislochar][2][i] - MatData[cislochar][2][i-1]))/(MatData[cislochar][0][i] - MatData[cislochar][0][i-1]);
01282 }
01283 }
01284 }
01285 kvyhl = 2;
01286 z2 = MatData[cislochar][1][prad];
01287 break;
01288 }
01289 case 30:{
01290 y = MatFunce[cislochar][0];
01291 z = MatFunce[cislochar][1];
01292 kvyhl = 30;
01293 z2 = -10000;
01294 break;
01295 }
01296 case 31:{
01297 y = MatFunce[cislochar][0];
01298 z = MatFunce[cislochar][1];
01299 if (cislochar == 6)
01300 {
01301 z2 = MatFunce[cislochar][2];
01302 }
01303 else {
01304 z2 = -10000;
01305 }
01306 kvyhl = 31;
01307 break;
01308 }
01309 case 32:{
01310 y = MatFunce[cislochar][0];
01311 z = MatFunce[cislochar][1];
01312 kvyhl = 32;
01313 z2 = -10000;
01314 break;
01315 }
01316 case 33:{
01317 kvyhl = 33;
01318
01319 z2 = -10000;
01320 z = -10000;
01321 y = -10000;
01322 break;
01323 }
01324 case 40:{
01325 int pradA, pradD;
01326 pradA = int(MatData[cislochar][0][0]);
01327 pradD = int(MatData[cislochar][2][0]);
01328 if (x <=0)
01329 {
01330 x = 0;
01331 y = MatData[cislochar][1][1];
01332 z = MatData[cislochar][3][1];
01333 }
01334 else
01335 {
01336 if (x > in)
01337 {
01338 y = MatData[cislochar][1][pradA];
01339 z = MatData[cislochar][3][pradD];
01340 }
01341 else
01342 {
01343 i = 1;
01344 while (x > MatData[cislochar][0][i]){
01345 i++;
01346 }
01347 y = MatData[cislochar][1][i-1] + ((x - MatData[cislochar][0][i-1])*(MatData[cislochar][1][i] - MatData[cislochar][1][i-1]))/(MatData[cislochar][0][i] - MatData[cislochar][0][i-1]);
01348 i = 1;
01349 while (x > MatData[cislochar][2][i]){
01350 i++;
01351 }
01352 z = MatData[cislochar][3][i-1] + ((x - MatData[cislochar][2][i-1])*(MatData[cislochar][3][i] - MatData[cislochar][3][i-1]))/(MatData[cislochar][2][i] - MatData[cislochar][2][i-1]);
01353 }
01354 }
01355 z2 = MatData[cislochar][1][pradA];
01356 kvyhl = 40;
01357 }
01358 break;
01359 }
01360 }
01361
01362 void kunmat2::sorption_izoterms_values(int kod, long ipp, double x1,double xpv, double ineq1, double & w, double & dwdf)
01363 {
01364 int s;
01365 double hmrh, mhmc, smc, u, a, n;
01366 double x2;
01367
01368 x2 = -10000.0;
01369
01370 if (x1 < 0) x1 = 0;
01371
01372
01373
01374 switch (kod){
01375 case 2:{
01376 if(x1>0.976)
01377 {
01378
01379 CorD(7,kd,0,x1,smc,a2,a3);
01380 CorD(6,kd,0,x1,a3,a2,mhmc);
01381 give_data(x1,mhmc, smc, 0.976,w);
01382 dwdf =1000 * (smc - mhmc)/(1 - 0.976);
01383 }
01384 else
01385 {
01386 CorD(6,s,1,x1,w,dwdf,a3);
01387 dwdf =derivation_sorption_izoterm_data(x1,x2,0,ipp)*1000;
01388 }
01389 break;
01390 }
01391 case 30:{
01392 CorD(7,s,0,x1,smc,a2,a3);
01393 CorD(6,s,0,x1,mhmc,hmrh,a3);
01394 give_data(x1,mhmc, smc, hmrh,w);
01395 if (x1 < hmrh)
01396 {
01397 dwdf = DerivativeOfTheSorptionIsotherm (x1, x2);
01398 }
01399 else {
01400 dwdf = DerivativeOfTheRetentionCurve(x1, x2);
01401 }
01402 break;
01403 }
01404 case 31:{
01405 CorD(6,kod,0,0,u,a,n);
01406 w = si_kk_hansen(x1, 0,u,a,n);
01407 dwdf = soptionizothermDerivation(x1, x2);
01408 break ;
01409 }
01410
01411 case 40:{
01412
01413 if(x1>=0.976)
01414 {
01415 CorD(7,kd,0,x1,smc,a2,a3);
01416 CorD(6,kd,0,x1,a1,a2,mhmc);
01417 give_data(x1,mhmc, smc, 0.976,w);
01418 dwdf =1000 * (smc - mhmc)/(1 - 0.976);
01419 }
01420 else{
01421 hystereze2 (6,x1, xpv, ineq1,w,dwdf, ipp);
01422 }
01423 break;
01424 }
01425
01426 default:{
01427 fprintf (stderr,"\n\n unknown definition of Sorption isotherm is required");
01428 fprintf (stderr,"\n in function sorption_izoterms_giva_data (file %s, line %d)\n",__FILE__,__LINE__);
01429 }
01430 }
01431
01432 }
01433
01434 double kunmat2::si_kk_hansen(double x1, double x2,double u, double a, double n)
01435 {
01436 return (u*pow((1-(log(x1))/a),(-1/n)));
01437 }
01438
01439 double kunmat2::derivation_sorption_izoterm_data(double x1, double x2, double x3, long ipp)
01440 {
01441 int i = 1;
01442 double derbi;
01443 if ( x1 < 0)
01444 {
01445 x1 = 0;
01446 }
01447
01448 while (x1 > MatData[6][0][i]){
01449 i++;
01450 }
01451 derbi = (MatData [6][1][i] - MatData [6][1][i-1])/(MatData [6][0][i] - MatData [6][0][i-1]);
01452 return(derbi);
01453 }
01454
01455
01456 void kunmat2::hystereze2 (int matchar, double x, double xpv, double ineq1,double & outvalue,double & outvalue2, long ipp)
01457 {
01458 double xav,yav,newder,ypv,yavA,yavD,ypvA,ypvD,derpvA,derpvD,other1,deravA,deravD;
01459 int hyst, other;
01460
01461 double t;
01462 t = Tp->time;
01463
01464
01465 if (t < Init[matchar])
01466 {
01467 switch (matchar)
01468 {
01469 case 5: {
01470 CorD(5,other,1,x,outvalue,outvalue2,other1);
01471
01472 }
01473 break;
01474 case 6: {
01475 CorD(6,other,1,x,outvalue,outvalue2,other1);
01476 outvalue2 = derivation_dy_dx(6,x,0,1)*1000;
01477 }
01478 break;
01479 }
01480 }
01481 else
01482 {
01483 switch(matchar)
01484 {
01485 case 5:{
01486 double yavANorm,yavDNorm,ypvANorm,ypvDNorm,derpvANorm,derpvDNorm,ypvNorm,yavNorm;
01487 double p, logP;
01488
01489 xav = x;
01490
01491 ypv = ineq1;
01492
01493 ypvNorm = log(ineq1);
01494 CorD(5,other,1,xav,yavA,yavD,other1);
01495 CorD(5,other,1,xpv,ypvA,ypvD,other1);
01496
01497 yavANorm = log(yavA);
01498
01499 yavDNorm = log(yavD);
01500
01501 ypvANorm = log(ypvA);
01502
01503 ypvDNorm = log(ypvD);
01504
01505 der_value_hyst(matchar,40,xpv,derpvA,derpvD,ipp);
01506
01507
01508 derpvANorm = log(derpvA);
01509 derpvDNorm = log(derpvD);
01510
01511
01512 if(xav > xpv)
01513 {
01514 hyst = 0;
01515 }
01516 else {
01517 hyst = 1;
01518 }
01519 switch(hyst)
01520 {
01521 case 0:{
01522 if(ypvNorm == ypvANorm)
01523 {
01524 yavNorm = yavANorm;
01525 }
01526 else {
01527 newder = (k1[matchar]*derpvDNorm*(ypvNorm-ypvANorm)*(ypvNorm-ypvANorm)+k2[matchar]*derpvANorm*(ypvNorm -ypvDNorm)*(ypvNorm - ypvDNorm))/((ypvDNorm-ypvANorm)*(ypvDNorm-ypvANorm));
01528
01529 p = -1*newder;
01530 logP = log(p);
01531 yavNorm =ypvNorm + logP*(xav-xpv)*k3[matchar];
01532
01533
01534 if(yavNorm > yavANorm) yavNorm = yavANorm;
01535 if(yavNorm < yavDNorm) yavNorm = yavDNorm;
01536 }
01537 break;
01538 }
01539 case 1:{
01540 if(ypvNorm == ypvDNorm)
01541 {
01542 yavNorm = yavDNorm;
01543 }
01544 else {
01545 newder = (k4[matchar]*derpvDNorm*(ypvNorm-ypvANorm)*(ypvNorm-ypvANorm)+k5[matchar]*derpvANorm*(ypvNorm -ypvDNorm)*(ypvNorm - ypvDNorm))/((ypvDNorm-ypvANorm)*(ypvDNorm-ypvANorm));
01546
01547
01548 p = -1*newder;
01549 logP = log(p);
01550 yavNorm =ypvNorm + logP*(xav-xpv)*k6[matchar];
01551
01552
01553 if(yavNorm < yavDNorm) yavNorm = yavDNorm;
01554 if(yavNorm > yavANorm) yavNorm = yavANorm;
01555 }
01556 break;
01557 }
01558 }
01559 yav = exp(yavNorm);
01560 outvalue = yav;
01561 outvalue2 = -10000.0;
01562 if (t>15559200)
01563 {
01564 if(t< 15724800)
01565 {
01566 if(ipp == 80)
01567 fprintf(Outt2,"%lf %e %e %e %e\n",t, xav,yavA,yav, yavD);
01568 }
01569 else{
01570 if(t> 31536000 && t< 31622400)
01571 {
01572 if(ipp == 80)
01573 fprintf(Outt2,"%lf %e %e %e %e\n",t, xav,yavA,yav, yavD);
01574 }
01575 }
01576
01577 }
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624 break;
01625 }
01626 case 6:
01627 {
01628 xav = x;
01629 ypv = ineq1;
01630 CorD(6,other,1,xav,yavA,yavD,other1);
01631 CorD(6,other,1,xpv,ypvA,ypvD,other1);
01632
01633
01634 der_value_hyst(matchar,40,xpv,derpvA,derpvD,ipp);
01635 if(xav >= xpv)
01636 {
01637 hyst = 0;
01638 }
01639 else{
01640 hyst = 1;
01641 }
01642
01643 switch(hyst)
01644 {
01645 case 0:{
01646 if(ypv == ypvA)
01647 {
01648 yav = yavA;
01649 }
01650 else {
01651 newder = (k1[matchar]*derpvD*(ypv-ypvA)*(ypv-ypvA)+k2[matchar]*derpvA*(ypv - ypvD)*(ypv - ypvD))/((ypvD-ypvA)*(ypvD-ypvA));
01652 yav = (xav-xpv)*newder*k3[matchar] + ypv;
01653 if(yav < yavA) yav = yavA;
01654 if(yav > yavD) yav = yavD;
01655 }
01656 break;
01657 }
01658 case 1:{
01659 if(ypv == ypvD)
01660 {
01661 yav = yavD;
01662 }
01663 else {
01664 newder = (k4[matchar]*derpvD*(ypv-ypvA)*(ypv-ypvA)+k5[matchar]*derpvA*(ypv -ypvD)*(ypv - ypvD))/((ypvD-ypvA)*(ypvD-ypvA));
01665 yav = (xav-xpv)*newder*k6[matchar] + ypv;
01666 if(yav < yavA) yav = yavA;
01667 if(yav > yavD) yav = yavD;
01668
01669 }
01670 break;
01671 }
01672 }
01673
01674 der_value_hyst(matchar,40,xav,deravA,deravD,ipp);
01675
01676 if (yav == yavA)
01677 newder = deravA;
01678 if (yav == yavD)
01679 newder = deravD;
01680
01681 double newder2;
01682
01683 newder2 = (yav-ypv)/(xav-xpv);
01684 outvalue = yav;
01685 outvalue2 = newder*1000;
01686 if (t>15559200)
01687 {
01688 if(t< 15724800)
01689 {
01690 if(ipp == 80){
01691
01692 }
01693 }
01694 else{
01695 if(t> 31536000 && t< 31622400)
01696 {
01697 if(ipp == 80){
01698
01699 }
01700 }
01701 }
01702
01703 }
01704
01705 }
01706 break;
01707 }
01708 }
01709 }
01710
01711 double kunmat2::derivation_dy_dx (int matchar, double prom, int pomk1, int pomk2)
01712
01713 {
01714 int i = 1;
01715 double outvalue;
01716
01717 while (prom > MatData[matchar][pomk1][i]){
01718 i++;
01719 }
01720 outvalue = (MatData [matchar][pomk2][i] - MatData [matchar][pomk2][i-1])/(MatData [matchar][pomk1][i] - MatData [matchar][pomk1][i-1]);
01721
01722 return(outvalue);
01723 }
01724
01725 void kunmat2::der_value_hyst (int matchar,int kod, double pv, double & outvalue,double & outvalue2, long ipp)
01726 {
01727
01728 switch (kod){
01729 case 0:{
01730 break;
01731 }
01732 case 1:{
01733 break;
01734 }
01735 case 2:{
01736 break;
01737 }
01738 case 30:{
01739 break;
01740 }
01741 case 31:{
01742 break;
01743 }
01744 case 40:{
01745
01746 outvalue = derivation_dy_dx(matchar,pv,0,1);
01747 outvalue2 = derivation_dy_dx(matchar,pv,2,3);
01748
01749 }
01750 break;
01751 }
01752 }
01753
01754 void kunmat2::kapa_values (int kod, long ipp,double x1, double xpv, double ineq1, double &kapa)
01755 {
01756 double kapa1, kapa2, kapa3, x1w;
01757 double smc, a, n, x2;
01758 int kd1;
01759 x2 = 0.0;
01760
01761 CorD(7,kd1,1,x1,smc,a,n);
01762
01763 if( x1 > smc) x1 = smc;
01764
01765 switch (kod){
01766 case 0: kapa = 0.0;
01767 break;
01768 case 1: CorD(5,kod,1,x1,kapa,kapa1,kapa2);
01769 break;
01770 case 2: CorD(5,kod,1,x1,kapa,kapa1,kapa2);
01771 break;
01772 case 30:{
01773 CorD(5,kod,1,x1,kapa1,kapa2,kapa3);
01774 kapa = kapa_exp(kapa1,kapa2,x1, x2, ipp);
01775 break;
01776 }
01777 case 40: hystereze2 (5,x1,xpv,ineq1,kapa,kapa1,ipp);
01778 break;
01779 }
01780 }
01781
01782
01783
01784
01785
01786
01787
01788 void kunmat2::aux_values (long ipp,double *in,double *inp, double *ineq,double *out)
01789 {
01790 double x1,x2,x1p, x2p, ypv1,ypv2,ypv3,ypv4,ypv5, ypv0;
01791
01792 double kapaA,w,dwdf;
01793 double t;
01794 t = Tp->time;
01795 x1 = in[0];
01796 x2 = in[1];
01797
01798 x1p = inp[0];
01799 x2p = inp[1];
01800
01801
01802 ypv0 =ineq[0];
01803 ypv1 =ineq[1];
01804 ypv2 =ineq[2];
01805 ypv3 =ineq[3];
01806 ypv4 =ineq[4];
01807 ypv5 =ineq[5];
01808
01809 sorption_izoterms_values(MatChar[6],ipp,x1,x1p,ypv0,w,dwdf);
01810
01811 kapa_values(MatChar[5],ipp,w,ypv0,ypv2, kapaA);
01812 out[3] = kapaA;
01813
01814 out[0] = w;
01815 out[1] = dwdf;
01816 out[2] = MatConst[7];
01817 out[4] = -10000.0;
01818 out[5] = -10000.0;
01819
01820 if (ipp == 1)
01821 fprintf(Outt,"t= %lf je w == %3.2e dwdrh je %3.2e a kapa je %3.2e \n",t,w,dwdf,kapaA);
01822
01823 }
01824
01825 void kunmat2::initvalues (long ipp,long ido)
01826 {
01827 double x1,x2;
01828
01829 x1 = Tm->ip[ipp].av[0];
01830 x2 = Tm->ip[ipp].av[1];
01831 }
01832
01833
01834
01835
01836
01837
01838 void kunmat2::give_values (long ipp,double *av, double *pv, double *eq)
01839 {
01840 av[0] = Tm->ip[ipp].av[0];
01841 av[1] = Tm->ip[ipp].av[1];
01842
01843 pv[0] = Tm->ip[ipp].pv[0];
01844 pv[1] = Tm->ip[ipp].pv[1];
01845
01846
01847 eq[0] = Tm->ip[ipp].eqother[0];
01848 eq[1] = Tm->ip[ipp].eqother[1];
01849 eq[2] = Tm->ip[ipp].eqother[2];
01850 eq[3] = Tm->ip[ipp].eqother[3];
01851 eq[4] = Tm->ip[ipp].eqother[4];
01852 eq[5] = Tm->ip[ipp].eqother[5];
01853
01854
01855 }
01856
01857
01858
01859
01860
01861
01862
01863 void kunmat2::save_values (long ipp,double *out)
01864 {
01865 Tm->ip[ipp].eqother[0]=out[0];
01866 Tm->ip[ipp].eqother[1]=out[1];
01867 Tm->ip[ipp].eqother[2]=out[2];
01868 Tm->ip[ipp].eqother[3]=out[3];
01869 Tm->ip[ipp].eqother[4]=out[4];
01870 Tm->ip[ipp].eqother[5]=out[5];
01871
01872 }