source: inundation/ga/storm_surge/examples/beach.py @ 653

"""Example of the inundationmodel.

A wave of water is washed up ontp a hypothetical beach.

To run:

python beach.py
"""

######################
# Module imports 
#

import sys
from os import sep, path
sys.path.append('..'+sep+'pyvolution')

from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
     Transmissive_boundary, Time_boundary, Constant_height, Weir

from pmesh2domain import pmesh_to_domain_instance
from config import data_dir

######################
# Domain
name = 'beach'
print 'Creating domain from %s.tsh' %name
domain = pmesh_to_domain_instance(name + '.tsh', Domain)

#HACK to remove internal boundary (until we get a solution)
for b in domain.boundary.keys():
    if domain.boundary[b] == 'internal':
        del domain.boundary[b]

domain = Domain(domain.coordinates, domain.triangles, domain.boundary)
print "Number of triangles = ", len(domain)

domain.store = True
domain.set_name(name)
print "Output being written to " + data_dir + sep + \
      domain.filename + "." + domain.format


def x_slope(x, y):
    from Numeric import choose, greater, ones
    cut = 65
    
    return choose( greater(x, cut), (-x/3, -cut/3*ones(x.shape) )) + cut/3 - 5


def inside(p,polygon):
    return True

def friction(x, y):
    """Read friction from polygon file and set accordingly
    """

    from Numeric import array
    x = array(x).astype('f')
    y = array(y).astype('f')    
    
    #Get polygon
    fid = open('beach.poly')
    lines = fid.readlines()
    fid.close()
    polygon = []
    for line in lines:
        fields = line.split(',')
        polygon.append( [float(fields[0]), float(fields[1])] )

    print polygon

    z = 0*x
    for i in range(len(x)):
        if inside( (x[i], y[i]), polygon):
            z[i] = 1
        else:
            z[i] = 0

friction([0],[1])    

domain.set_quantity('elevation', x_slope)
domain.set_quantity('level', x_slope)
domain.set_quantity('friction', 0.07)

print domain.get_extent()
import sys; sys.exit()

######################
# Boundary conditions

print 'Boundaries'
Br = Reflective_boundary(domain)
Bt = Transmissive_boundary(domain)

#Constant inflow
Bd = Dirichlet_boundary([0.0, -1.0, 0.0])

print 'Available boundary tags are', domain.get_boundary_tags()

#Set boundary conditions
tags = {}
tags['ocean'] = Bd
tags['wall'] = Br 
tags['exterior'] = Br
tags['external'] = Br 
domain.set_boundary(tags)


domain.check_integrity()

######################
#Evolution
for t in domain.evolve(yieldstep = 0.5, finaltime = 100):
    domain.write_time()
    
print 'Done'