00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "richardscontam.h"
00013 #include "richardscontam.h"
00014 #include "globalt.h"
00015
00016
00017 richardscontam::richardscontam()
00018 {
00019
00020 alpha=0.0;
00021 n=0.0;
00022 m=0.0;
00023
00024 dim=0;
00025
00026
00027 kksxx=0.0;
00028 kksxy=0.0;
00029 kksxz=0.0;
00030 kksyy=0.0;
00031 kksyz=0.0;
00032 kkszz=0.0;
00033
00034 thetas=0.0;
00035
00036 thetar=0.0;
00037
00038 }
00039
00040 richardscontam::~richardscontam()
00041 {
00042
00043 }
00044
00045
00046
00047
00048
00049
00050
00051
00052 void richardscontam::read(XFILE *in)
00053 {
00054
00055
00056 dim = Tt->give_dimension (0);
00057 if (dim<1 || dim>3){
00058 print_err("wrong dimension of problem solved (%ld)",__FILE__,__LINE__,__func__,dim);
00059 }
00060
00061
00062
00063
00064
00065
00066 rich.read (in);
00067
00068 if (dim==1){
00069 xfscanf (in,"%lf",&kksxx);
00070 }
00071 if (dim==2){
00072 xfscanf (in,"%lf %lf %lf",&kksxx,&kksxy,&kksyy);
00073 }
00074 if (dim==3){
00075 xfscanf (in,"%lf %lf %lf %lf %lf %lf",&kksxx,&kksxy,&kksxz,&kksyy,&kksyz,&kkszz);
00076 }
00077
00078
00079 xfscanf (in,"%lf %lf %lf %lf %lf %lf",&alpha,&n,&m,&thetas,&thetar,&storage);
00080 }
00081
00082
00083
00084
00085
00086
00087
00088
00089 void richardscontam::print(FILE *out)
00090 {
00091 rich.print (out);
00092
00093 if (dim==1){
00094 fprintf (out," %lf\n",kksxx);
00095 }
00096 if (dim==2){
00097 fprintf (out," %lf %lf %lf\n",kksxx,kksxy,kksyy);
00098 }
00099 if (dim==3){
00100 fprintf (out," %lf %lf %lf %lf %lf %lf\n",kksxx,kksxy,kksxz,kksyy,kksyz,kkszz);
00101 }
00102
00103
00104
00105
00106
00107 fprintf (out," %lf %lf %lf %lf %lf %lf\n",alpha,n,m,thetas,thetar,storage);
00108 }
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121 void richardscontam::matcond (matrix &d,long ri,long ci,long ipp)
00122 {
00123 long n;
00124 n = d.n;
00125
00126 switch (n){
00127 case 1:{
00128 matcond1d (d,ri,ci,ipp);
00129 break;
00130 }
00131 case 2:{
00132 matcond2d (d,ri,ci,ipp);
00133 break;
00134 }
00135 case 3:{
00136 matcond3d (d,ri,ci,ipp);
00137 break;
00138 }
00139 default:{
00140 print_err("unknown number of flux components is required",__FILE__,__LINE__,__func__);
00141 }
00142 }
00143 }
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 void richardscontam::matcond1d (matrix &d,long ri,long ci,long ipp)
00155 {
00156 double h,k;
00157
00158 if (ri==0 && ci==0){
00159
00160
00161 h = Tm->ip[ipp].av[0];
00162
00163 k = rich.kk_value (h);
00164 }
00165 if (ri==0 && ci==1){
00166
00167 k=0.0;
00168 }
00169 if (ri==1 && ci ==0){
00170
00171 k=0.0;
00172 }
00173 if (ri==1 && ci==1){
00174
00175
00176 h = Tm->ip[ipp].av[0];
00177
00178
00179 k = cc_value (h);
00180
00181 }
00182
00183
00184 d[0][0] = k*kksxx;
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 void richardscontam::matcond2d (matrix &d,long ri,long ci,long ipp)
00196 {
00197 double h,k;
00198
00199 if (ri==0 && ci==0){
00200
00201
00202 h = Tm->ip[ipp].av[0];
00203
00204 k = rich.kk_value (h);
00205 }
00206 if (ri==0 && ci==1){
00207
00208 k=0.0;
00209 }
00210 if (ri==1 && ci ==0){
00211
00212 k=0.0;
00213 }
00214 if (ri==1 && ci==1){
00215
00216
00217 h = Tm->ip[ipp].av[0];
00218
00219
00220 k = cc_value (h);
00221
00222 }
00223
00224 d[0][0] = k*kksxx; d[0][1] = k*kksxy;
00225 d[1][0] = k*kksxy; d[1][1] = k*kksyy;
00226
00227 }
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 void richardscontam::matcond3d (matrix &d,long ri,long ci,long ipp)
00239 {
00240 double h,k;
00241
00242 if (ri==0 && ci==0){
00243
00244
00245 h = Tm->ip[ipp].av[0];
00246
00247 k = rich.kk_value (h);
00248 }
00249 if (ri==0 && ci==1){
00250
00251 k=0.0;
00252 }
00253 if (ri==1 && ci ==0){
00254
00255 k=0.0;
00256 }
00257 if (ri==1 && ci==1){
00258
00259
00260 h = Tm->ip[ipp].av[0];
00261
00262
00263 k = cc_value (h);
00264
00265 }
00266
00267 d[0][0] = k*kksxx; d[0][1] = k*kksxy; d[0][2]=k*kksxz;
00268 d[1][0] = k*kksxy; d[1][1] = k*kksyy; d[1][2]=k*kksyz;
00269 d[2][0] = k*kksxz; d[2][1] = k*kksyz; d[2][2]=k*kkszz;
00270 }
00271
00272 void richardscontam::matcond2 (matrix &d,long ri,long ci,long ipp)
00273 {
00274 long n;
00275 n = d.n;
00276
00277 switch (n){
00278 case 1:{
00279 matcond1d_2 (d,ri,ci,ipp);
00280 break;
00281 }
00282 case 2:{
00283 matcond2d_2 (d,ri,ci,ipp);
00284 break;
00285 }
00286 case 3:{
00287 matcond3d_2 (d,ri,ci,ipp);
00288 break;
00289 }
00290 default:{
00291 print_err("unknown number of components is required",__FILE__,__LINE__,__func__);
00292 }
00293 }
00294 }
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 void richardscontam::matcond1d_2 (matrix &d,long ri,long ci,long ipp)
00307 {
00308 double h;
00309
00310 fillm(0.0,d);
00311 if (ri==0 && ci==0){
00312
00313 h = Tm->ip[ipp].av[0];
00314
00315 d[0][0]=0.0-rich.dkkdh_value (h);
00316 }
00317 }
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331 void richardscontam::matcond2d_2 (matrix &d,long ri,long ci,long ipp)
00332 {
00333 double h;
00334
00335 fillm(0.0,d);
00336 if (ri==0 && ci==0){
00337
00338 h = Tm->ip[ipp].av[0];
00339
00340 d[0][0]=0.0;
00341 d[0][1]=0.0-rich.dkkdh_value (h);
00342 }
00343 }
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 void richardscontam::matcond3d_2 (matrix &d,long ri,long ci,long ipp)
00356 {
00357 double h;
00358
00359 fillm(0.0,d);
00360
00361
00362 h = Tm->ip[ipp].av[0];
00363 if (ri==0 && ci==0){
00364 d[0][0]=0.0;
00365 d[0][1]=0.0;
00366 d[0][2]=0.0-rich.dkkdh_value (h);
00367 }
00368 }
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380 void richardscontam::matcap (double &c,long ri,long ci,long ipp)
00381 {
00382 double h;
00383
00384
00385 h = Tm->ip[ipp].av[0];
00386
00387
00388
00389
00390 }
00391
00392 double richardscontam::cc_value (double h)
00393 {
00394 return 0.0;
00395 }