00001
00002
00003
00004
00005
00006
00007
00008
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <math.h>
00012 #include "constrel.h"
00013 #include "twomediac.h"
00014 #include "coupmatu.h"
00015 #include "coupmatl.h"
00016 #include "globalc.h"
00017 #include "global.h"
00018 #include "globalt.h"
00019 #include "consol_awf2c.h"
00020 #include "intpoints.h"
00021
00022
00023 con_awf2matc::con_awf2matc()
00024 {
00025 alpha = 1.0;
00026
00027 rhos = 2000.0;
00028
00029 rhow = 1000.0;
00030
00031 phi0 = 0.297;
00032
00033 }
00034
00035 con_awf2matc::~con_awf2matc()
00036 {}
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 void con_awf2matc::read(XFILE *in)
00047 {
00048 xfscanf (in,"%lf %lf %lf",&rhos, &alpha, &phi0);
00049 }
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 void con_awf2matc::print(FILE *out)
00060 {
00061 fprintf (out," %lf %lf %lf\n",rhos,alpha,phi0);
00062 }
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 void con_awf2matc::matcond_u (matrix &d,long ri,long ci,long ipp)
00076 {
00077 long m;
00078 m = d.m;
00079
00080 switch (m){
00081 case 1:{
00082 matcond1d_u (d,ri,ci,ipp);
00083 break;
00084 }
00085 case 2:{
00086 matcond2d_u (d,ri,ci,ipp);
00087 break;
00088 }
00089 case 3:{
00090 matcond3d_u (d,ri,ci,ipp);
00091 break;
00092 }
00093 default:{
00094 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00095 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00096 }
00097 }
00098 }
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 void con_awf2matc::matcond1d_u (matrix &d,long ri,long ci,long ipp)
00113 {
00114 double k;
00115 double pw,pg;
00116 k = 0.0;
00117
00118 pw = Cmu->ip[ipp].av[0];
00119 pg = Cmu->ip[ipp].av[1];
00120
00121 if((ri == 0) && (ci == 0))
00122 k = get_kuw(pw,pg);
00123 if((ri == 0) && (ci == 1))
00124 k = get_kug(pw,pg);
00125
00126 d[0][0] = k;
00127 }
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 void con_awf2matc::matcond2d_u (matrix &d,long ri,long ci,long ipp)
00141 {
00142 double k;
00143 double pw,pg;
00144 k = 0.0;
00145
00146 pw = Cmu->ip[ipp].av[0];
00147 pg = Cmu->ip[ipp].av[1];
00148
00149 if((ri == 0) && (ci == 0))
00150 k = get_kuw(pw,pg);
00151 if((ri == 0) && (ci == 1))
00152 k = get_kug(pw,pg);
00153
00154 fillm(0.0,d);
00155
00156 d[0][0] = k; d[0][1] = 0.0;
00157 d[1][0] = 0.0; d[1][1] = k;
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 void con_awf2matc::matcond3d_u (matrix &d,long ri,long ci,long ipp)
00172 {
00173 double k;
00174 double pw,pg;
00175 k = 0.0;
00176
00177 pw = Cmu->ip[ipp].av[0];
00178 pg = Cmu->ip[ipp].av[1];
00179
00180 if((ri == 0) && (ci == 0))
00181 k = get_kuw(pw,pg);
00182 if((ri == 0) && (ci == 1))
00183 k = get_kug(pw,pg);
00184
00185 fillm(0.0,d);
00186
00187 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00188 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00189 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204 void con_awf2matc::matcap_u (matrix &d,long ri,long ci,long ipp)
00205 {
00206 long m;
00207
00208 m = d.m;
00209
00210 switch (m){
00211 case 1:{
00212 matcap1d_u (d,ri,ci,ipp);
00213 break;
00214 }
00215 case 3:{
00216 matcap2d_u (d,ri,ci,ipp);
00217 break;
00218 }
00219 case 6:{
00220 matcap3d_u (d,ri,ci,ipp);
00221 break;
00222 }
00223 default:{
00224 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00225 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00226 }
00227 }
00228 }
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 void con_awf2matc::matcap1d_u (matrix &d,long ri,long ci,long ipp)
00243 {
00244 double pw,pg,c;
00245
00246 pw = Cmu->ip[ipp].av[0];
00247 pg = Cmu->ip[ipp].av[1];
00248
00249 if((ri == 0) && (ci == 0))
00250 c = get_capuw(pw,pg);
00251 if((ri == 0) && (ci == 1))
00252 c = get_capug(pw,pg);
00253
00254 d[0][0] = c;
00255 }
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 void con_awf2matc::matcap2d_u (matrix &d,long ri,long ci,long ipp)
00270 {
00271 double pw,pg,c;
00272
00273 pw = Cmu->ip[ipp].av[0];
00274 pg = Cmu->ip[ipp].av[1];
00275
00276 if((ri == 0) && (ci == 0))
00277 c = get_capuw(pw,pg);
00278 if((ri == 0) && (ci == 1))
00279 c = get_capug(pw,pg);
00280
00281
00282 fillm(0.0,d);
00283
00284 d[0][0] = c;
00285 d[1][0] = c;
00286 d[2][0] = 0.0;
00287 }
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301 void con_awf2matc::matcap3d_u (matrix &d,long ri,long ci,long ipp)
00302 {
00303 double pw,pg,c;
00304
00305 pw = Cmu->ip[ipp].av[0];
00306 pg = Cmu->ip[ipp].av[1];
00307
00308 if((ri == 0) && (ci == 0))
00309 c = get_capuw(pw,pg);
00310 if((ri == 0) && (ci == 1))
00311 c = get_capug(pw,pg);
00312
00313
00314 fillm(0.0,d);
00315
00316 d[0][0] = c;
00317 d[1][0] = c;
00318 d[2][0] = c;
00319 d[3][0] = 0.0;
00320 d[4][0] = 0.0;
00321 d[5][0] = 0.0;
00322 }
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336 void con_awf2matc::matcond_l (matrix &d,long ri,long ci,long ipp)
00337 {
00338 long m;
00339 m = d.n;
00340
00341 switch (m){
00342 case 1:{
00343 matcond1d_l (d,ri,ci,ipp);
00344 break;
00345 }
00346 case 2:{
00347 matcond2d_l (d,ri,ci,ipp);
00348 break;
00349 }
00350 case 3:{
00351 matcond3d_l (d,ri,ci,ipp);
00352 break;
00353 }
00354 default:{
00355 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00356 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00357 }
00358 }
00359 }
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 void con_awf2matc::matcond1d_l (matrix &d,long ri,long ci,long ipp)
00374 {
00375 double k;
00376 double pw,pg;
00377 k = 0.0;
00378
00379 pw = Cml->ip[ipp].av[0];
00380 pg = Cml->ip[ipp].av[1];
00381
00382 if((ri == 0) && (ci == 0))
00383 k = get_kwu(pw,pg);
00384 if((ri == 0) && (ci == 1))
00385 k = get_kgu(pw,pg);
00386
00387 d[0][0] = k;
00388 }
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 void con_awf2matc::matcond2d_l (matrix &d,long ri,long ci,long ipp)
00402 {
00403 double k;
00404 double pw,pg;
00405 k = 0.0;
00406
00407 pw = Cml->ip[ipp].av[0];
00408 pg = Cml->ip[ipp].av[1];
00409
00410 if((ri == 0) && (ci == 0))
00411 k = get_kwu(pw,pg);
00412 if((ri == 0) && (ci == 1))
00413 k = get_kgu(pw,pg);
00414
00415 fillm(0.0,d);
00416
00417 d[0][0] = k; d[0][1] = 0.0;
00418 d[1][0] = 0.0; d[1][1] = k;
00419 }
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 void con_awf2matc::matcond3d_l (matrix &d,long ri,long ci,long ipp)
00433 {
00434 double k;
00435 double pw,pg;
00436 k = 0.0;
00437
00438 pw = Cml->ip[ipp].av[0];
00439 pg = Cml->ip[ipp].av[1];
00440
00441 if((ri == 0) && (ci == 0))
00442 k = get_kwu(pw,pg);
00443 if((ri == 0) && (ci == 1))
00444 k = get_kgu(pw,pg);
00445
00446 fillm(0.0,d);
00447
00448 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00449 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00450 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00451 }
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 void con_awf2matc::matcap_l (matrix &d,long ri,long ci,long ipp)
00466 {
00467 long m;
00468
00469 m = d.n;
00470
00471 switch (m){
00472 case 1:{
00473 matcap1d_l (d,ri,ci,ipp);
00474 break;
00475 }
00476 case 3:{
00477 matcap2d_l (d,ri,ci,ipp);
00478 break;
00479 }
00480 case 6:{
00481 matcap3d_l (d,ri,ci,ipp);
00482 break;
00483 }
00484 default:{
00485 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00486 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00487 }
00488 }
00489 }
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503 void con_awf2matc::matcap1d_l (matrix &d,long ri,long ci,long ipp)
00504 {
00505 double pw,pg,c;
00506 c = 0.0;
00507
00508 pw = Cml->ip[ipp].av[0];
00509 pg = Cml->ip[ipp].av[1];
00510
00511 if((ri == 0) && (ci == 0))
00512 c = get_capwu(pw,pg);
00513 if((ri == 0) && (ci == 1))
00514 c = get_capgu(pw,pg);
00515
00516 d[0][0] = c;
00517 }
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531 void con_awf2matc::matcap2d_l (matrix &d,long ri,long ci,long ipp)
00532 {
00533 double pw,pg,c;
00534 c = 0.0;
00535
00536 pw = Cml->ip[ipp].av[0];
00537 pg = Cml->ip[ipp].av[1];
00538
00539 if((ri == 0) && (ci == 0))
00540 c = get_capwu(pw,pg);
00541 if((ri == 0) && (ci == 1))
00542 c = get_capgu(pw,pg);
00543
00544
00545 fillm(0.0,d);
00546
00547 d[0][0] = c;
00548 d[1][0] = c;
00549 d[2][0] = 0.0;
00550 }
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564 void con_awf2matc::matcap3d_l (matrix &d,long ri,long ci,long ipp)
00565 {
00566 double pw,pg,c;
00567 c = 0.0;
00568
00569 pw = Cml->ip[ipp].av[0];
00570 pg = Cml->ip[ipp].av[1];
00571
00572 if((ri == 0) && (ci == 0))
00573 c = get_capwu(pw,pg);
00574 if((ri == 0) && (ci == 1))
00575 c = get_capgu(pw,pg);
00576
00577
00578 fillm(0.0,d);
00579
00580 d[0][0] = c;
00581 d[1][0] = c;
00582 d[2][0] = c;
00583 d[3][0] = 0.0;
00584 d[4][0] = 0.0;
00585 d[5][0] = 0.0;
00586 }
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599 void con_awf2matc::rhs_volume (matrix &d,long ri,long ci,long ipp)
00600 {
00601 long m;
00602 m = d.m;
00603
00604 switch (m){
00605 case 1:{
00606 rhs1d1 (d,ri,ci,ipp);
00607 break;
00608 }
00609 case 2:{
00610 rhs2d1 (d,ri,ci,ipp);
00611 break;
00612 }
00613 case 3:{
00614 rhs3d1 (d,ri,ci,ipp);
00615 break;
00616 }
00617 default:{
00618 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00619 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00620 }
00621 }
00622 }
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633 void con_awf2matc::rhs1d1 (matrix &d,long ri,long ci,long ipp)
00634 {
00635 double f,pw,pg;
00636 double g;
00637
00638 pw = Cmu->ip[ipp].av[0];
00639 pg = Cmu->ip[ipp].av[1];
00640
00641 f = get_fu1(pw,pg);
00642 g = Tp->gr1;
00643
00644 fillm(0.0,d);
00645 d[0][0] = f*g;
00646 }
00647
00648
00649
00650
00651
00652
00653
00654
00655 void con_awf2matc::rhs2d1 (matrix &d,long ri,long ci,long ipp)
00656 {
00657 double f,pw,pg;
00658 double *g;
00659 g = new double [2];
00660
00661 pw = Cmu->ip[ipp].av[0];
00662 pg = Cmu->ip[ipp].av[1];
00663
00664 f = get_fu1(pw,pg);
00665
00666 g[0] = Tp->gr1;
00667 g[1] = Tp->gr2;
00668
00669 fillm(0.0,d);
00670 d[0][0] = f*g[0];
00671 d[1][0] = f*g[1];
00672 }
00673
00674
00675
00676
00677
00678
00679
00680
00681 void con_awf2matc::rhs3d1 (matrix &d,long ri,long ci,long ipp)
00682 {
00683 double f,pw,pg;
00684 double *g;
00685 g = new double [3];
00686
00687 pw = Cmu->ip[ipp].av[0];
00688 pg = Cmu->ip[ipp].av[1];
00689
00690 f = get_fu1(pw,pg);
00691
00692 g[0] = Tp->gr1;
00693 g[1] = Tp->gr2;
00694 g[2] = Tp->gr3;
00695
00696 fillm(0.0,d);
00697 d[0][0] = f*g[0];
00698 d[1][0] = f*g[1];
00699 d[2][0] = f*g[2];
00700 }
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712 double con_awf2matc::get_sw(double pw, double pg)
00713 {
00714 double sw,pc;
00715
00716
00717 pc = -pw;
00718
00719 sw = 1.0 - (1.9722e-11)*pow(pc,2.4279);
00720
00721 return(sw);
00722 }
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734 double con_awf2matc::get_kuw(double pw, double pg)
00735 {
00736 double sw,kuw;
00737
00738 sw = get_sw(pw,pg);
00739
00740 kuw = -sw*alpha;
00741
00742 return(kuw);
00743
00744 }
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757 double con_awf2matc::get_kug(double pw, double pg)
00758 {
00759 double sw,sg,kug;
00760
00761 sw = get_sw(pw,pg);
00762 sg = 1.0 - sw;
00763
00764 kug = -sg*alpha;
00765
00766 kug = 0.0;
00767
00768 return(kug);
00769
00770 }
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783 double con_awf2matc::get_kgu(double pw, double pg)
00784 {
00785 double kgu;
00786
00787 kgu = 0.0;
00788
00789 return(kgu);
00790
00791 }
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803 double con_awf2matc::get_kwu(double pw, double pg)
00804 {
00805 double kwu;
00806
00807 kwu = 0.0;
00808
00809 return(kwu);
00810
00811 }
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823 double con_awf2matc::get_capuw(double pw, double pg)
00824 {
00825 double capuw;
00826
00827 capuw = 0.0;
00828
00829 return(capuw);
00830
00831 }
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843 double con_awf2matc::get_capwu(double pw, double pg)
00844 {
00845 double sw,capwu;
00846
00847 sw = get_sw(pw,pg);
00848
00849 capwu = sw*alpha;
00850
00851 return(capwu);
00852
00853 }
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865 double con_awf2matc::get_capug(double pw, double pg)
00866 {
00867 double capug;
00868
00869 capug = 0.0;
00870
00871 return(capug);
00872
00873 }
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886 double con_awf2matc::get_capgu(double pw, double pg)
00887 {
00888 double sg,sw,capgu;
00889
00890 sw = get_sw(pw,pg);
00891 sg = 1.0-sw;
00892
00893 capgu = sg*alpha;
00894
00895 return(capgu);
00896
00897 }
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910 double con_awf2matc::get_fu1(double pw, double pg)
00911 {
00912 double rhog,sw,sg,n,fu1;
00913
00914 sw = get_sw(pw,pg);
00915 sg = 1.0 - sw;
00916 n = phi0;
00917 rhog = 1.25;
00918
00919 fu1 = rhos*(n-1) + sw*n*rhow + sg*n*rhog;
00920
00921
00922 fu1 = 0.0;
00923
00924 return(fu1);
00925
00926 }