source: anuga_work/development/dam_test_from_brad_2007/dam_sample.py @ 4631

"""
Dam break!
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------

from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross

from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary

from anuga.utilities.polygon import read_polygon, plot_polygons
from anuga.utilities.polygon import Polygon_function

#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
length = 20.
width = 10.
dx = dy = 2           # Resolution: Length of subdivisions on both axes

points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
                                               len1=length, len2=width)

domain = Domain(points, vertices, boundary)   # Create domain
domain.set_name('dam_sample')                 # Output name
#domain.set_store_vertices_uniquely(True)      # Look at the 'real' triangles
domain.beta_h = 0.0
domain.tight_slope_limiters = 1
domain.set_minimum_allowed_height(0.01)

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------

# Polygons
study = [[0,0], [20,0], [20,10], [0,10], [0,0]]
#dam = [[0,0], [6,0], [6,10], [0,10], [0,0]]
poly_dam = read_polygon('dam.csv')
print poly_dam
#plot_polygons([poly_dam])


def dam1(x,y):
    z = -x/10                 #initialise with slope (this is required to create the vector z?)
    X = len(x)                  #find the number of vector elements
    for i in range(X):
        if 6<= x[i]<=7 or x[i]<1:    #dam wall
                z[i] = 2
        if 6<= x[i] <=7 and 4<= y[i] <=6: #remove unwanted wall
                z[i] = -x[i]/10
    return z

def water(x,y):
    z = -x/10                 #initialise with slope (this is required to create the vector z)
    X = len(x)                #find the number of vector elements

    for i in range(X):
        if x[i]<=6:    #dam
            z[i]= 2
            
    return z

##manual definition
#
#domain.set_quantity('elevation', water)                # Use function for elevation
#domain.set_quantity('friction', 0.01)                       # Constant friction 
#domain.set_quantity('stage', expression='elevation')        # take elevation
#domain.set_quantity('elevation', dam1)                # Use function for elevation


# Assign values to quantities
domain.set_quantity('elevation', dam1)                # Use function for elevation
domain.set_quantity('friction', 0.01)                 # Constant friction 
domain.set_quantity('stage', Polygon_function( [(study, -1),(poly_dam, 2)] )) 


#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
Bi = Dirichlet_boundary([1.5, 0, 0])         # Inflow
Bo = Dirichlet_boundary([-20/10, 0, 0])      # Outflow
Br = Reflective_boundary(domain)             # Solid reflective wall

domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})


#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 0.2, finaltime = 20):
    domain.write_time(track_speeds=False)