Changeset 7835


Ignore:
Timestamp:
Jun 14, 2010, 9:50:26 PM (15 years ago)
Author:
steve
Message:
 
Location:
anuga_work/development/flow_1d
Files:
1 added
1 deleted
1 moved

Legend:

Unmodified
Added
Removed

  • anuga_work/development/flow_1d/generic_1d/util.py

    r7834 r7835  
    5656        return 0.0
    5757
    58 def calculate_wetted_area(x1,x2,z1,z2,w1,w2):
    59     if (w1 > z1) & (w2 < z2) & (z1 <= z2):
    60         x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
    61         A = 0.5*(w1-z1)*(x-x1)
    62         L = x-x1
    63     elif (w1 < z1) & (w2 > z2) & (z1 < z2):
    64         x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
    65         A = 0.5*(w2-z2)*(x2-x)
    66         L = x2-x
    67     elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
    68         x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
    69         A = 0.5*(w2-z2)*(x2-x)
    70         L = x2-x
    71     elif (w1 > z1) & (w2 < z2) & (z1 > z2):
    72         x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
    73         A = 0.5*(w1-z1)*(x-x1)
    74         L = x-x1
    75     elif (w1 <= z1) & (w2 <= z2):
    76         A = 0.0
    77     elif (w1 == z1) & (w2 > z2) & (z2 < z1):
    78         A = 0.5*(x2-x1)*(w2-z2)
    79     elif (w2 == z2) & (w1 > z1) & (z1 < z2):
    80         A = 0.5*(x2-x1)*(w1-z1)
    81     return A
    82 
    83 
    84 def calculate_new_wet_area(x1,x2,z1,z2,A):
    85     from numpy import sqrt
    86     min_centroid_height = 1.0e-3
    87     # Assumes reconstructed stage flat in a wetted cell
    88     M = (z2-z1)/(x2-x1)
    89     L = (x2-x1)
    90     min_area = min_centroid_height*L
    91     max_area = 0.5*(x2-x1)*abs(z2-z1)
    92     if A < max_area:
    93         if (z1 < z2):
    94             x = sqrt(2*A/M)+x1
    95             wet_len = x-x1
    96             wc = z1 + sqrt(M*2*A)
    97         elif (z2 < z1):
    98             x = -sqrt(-2*A/M)+x2
    99             wet_len = x2-x
    100             wc = z2+sqrt(-M*2*A)
    101         else:
    102             wc = A/L+0.5*(z1+z2)
    103             wet_len = x2-x1
    104     else:
    105         wc = 0.5*(z1+z2)+A/L
    106         wet_len = x2-x1
    107            
    108     return wc,wet_len
    109 
    110 def calculate_new_wet_area_analytic(x1,x2,z1,z2,A,t):
    111     min_centroid_height = 1.0e-3
    112     # Assumes reconstructed stage flat in a wetted cell
    113     M = (z2-z1)/(x2-x1)
    114     L = (x2-x1)
    115     min_area = min_centroid_height*L
    116     max_area = 0.5*(x2-x1)*abs(z2-z1)
    117     w1,uh1 = analytic_cannal(x1,t)
    118     w2,uh2 = analytic_cannal(x2,t)
    119     if (w1 > z1) & (w2 < z2) & (z1 <= z2):
    120         print "test1"
    121         x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
    122         wet_len = x-x1
    123     elif (w1 < z1) & (w2 > z2) & (z1 < z2):
    124         print "test2"
    125         x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
    126         wet_len = x2-x
    127     elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
    128         print "test3"
    129         x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
    130         wet_len = x2-x
    131     elif (w1 > z1) & (w2 < z2) & (z1 > z2):
    132         print "test4"
    133         x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
    134         wet_len = x-x1
    135     elif (w1 >= z1) & (w2 >= z2):
    136         print "test5"
    137         wet_len = x2-x1
    138     else: #(w1 <= z1) & (w2 <= z2)
    139         print "test5"
    140         if (w1 > z1) | (w2 > z2):
    141             print "ERROR"
    142         wet_len = x2-x1       
    143     return w1,w2,wet_len,uh1,uh2
    144 
    145 def analytic_cannal(C,t):
    146 
    147     import math
    148     #N = len(C)
    149     #u = zeros(N,numpy.float)    ## water velocity
    150     #h = zeros(N,numpy.float)    ## water depth
    151     x = C
    152     g = 9.81
    153 
    154 
    155     ## Define Basin Bathymetry
    156     #z_b = zeros(N,numpy.float) ## elevation of basin
    157     #z = zeros(N,numpy.float)   ## elevation of water surface
    158     z_infty = 10.0       ## max equilibrium water depth at lowest point.
    159     L_x = 2500.0         ## width of channel
    160 
    161     A0 = 0.5*L_x                  ## determines amplitudes of oscillations
    162     omega = math.sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
    163 
    164     x1 = A0*cos(omega*t)-L_x # left shoreline
    165     x2 = A0*cos(omega*t)+L_x # right shoreline
    166     if (x >=x1) & (x <= x2):
    167         z_b = z_infty*(x**2/L_x**2) ## or A0*cos(omega*t)\pmL_x
    168         u = -A0*omega*sin(omega*t)
    169         z = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t))
    170     else:
    171        z_b = z_infty*(x**2/L_x**2)
    172        u=0.0
    173        z = z_b
    174     h = z-z_b
    175     return z,u*h
Note: See TracChangeset for help on using the changeset viewer.